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Preface 


This quarterly publication provides archival reports on developments in programs 
managed by JPL’s Telecommunications and Mission Operations Directorate (TMOD), 
which now includes the former Telecommunications and Data Acquisition (TDA) Office. 
In space communications, radio navigation, radio science, and ground-based radio and 
radar astronomy, it reports on activities of the Deep Space Network (DSN) in planning, 
supporting research and technology, implementation, and operations. Also included are 
standards activity at JPL for space data and information systems and reimbursable 
DSN work performed for other space agencies through NASA. The preceding work is 
all performed for NASA’s Office of Space Communications (OSC). 

TMOD also performs work funded by other NASA program offices through and 
with the cooperation of OSC. The first of these is the Orbital Debris Radar Program 
funded by the Office of Space Systems Development. It exists at Goldstone only and 
makes use of the planetary radar capability when the antennas are configured as science 
instruments making direct observations of the planets, their satellites, and asteroids of 
our solar system. The Office of Space Sciences funds the data reduction and science 
analyses of data obtained by the Goldstone Solar System Radar. The antennas at all 
three complexes are also configured for radio astronomy research and, as such, conduct 
experiments funded by the National Science Foundation in the U.S. and other agencies 
at the overseas complexes. These experiments are either in microwave spectroscopy or 
very long baseline interferometry. 

Finally, tasks funded under the JPL Director’s Discretionary Fund and the Caltech 
President’s Fund that involve TMOD are included. 

This and each succeeding issue of The Telecommunications and Data Acquisition 
Progress Report will present material in some, but not necessarily all, of the aforemen- 
tioned programs. 
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Orbit Determination of Highly Elliptical Earth 
Orbiters Using Improved Doppler 
Data-Processing Modes 

J. A. Estefan 
Navigation Systems Section 


A navigation error covariance analysis of four highly elliptical Earth orbits is de- 
scribed, with apogee heights ranging from 20,000 to 76,800 km and perigee heights 
ranging from 1,000 to 5,000 km. This analysis differs from earlier studies in that 
improved navigation data-processing modes were used to reduce the radio metric 
data. For this study, X-band (8.4-GHz) Doppler data were assumed to be acquired 
from two Deep Space Network radio antennas and reconstructed orbit errors propa- 
gated over a single day. Doppler measurements were formulated as total-count phase 
measurements and compared to the traditional formulation of differenced-count fre- 
quency measurements. In addition, an enhanced data-filtering strategy' was used, 
which treated the principal ground system calibration errors affecting the data as 
filter parameters. Results suggest that a 40- to 60-percent accuracy improvement 
may be achievable over traditional data-processing modes in reconstructed orbit 
errors , with a substantial reduction in reconstructed velocity errors at perigee. His- 
torically, this has been a regime in which stringent navigation requirements have 
been difficult to meet by conventional methods. 


I. Introduction 

The principal focus of recent navigation-related research has been on the development of new or im- 
proved navigation techniques to simultaneously improve performance, while reducing navigation-related 
requirements levied upon spacecraft and associated mission operations. The motivation for such a fo- 
cused effort is clear— tighter budgetary constraints imposed on current and future NASA space science 
missions. Advanced studies of interplanetary mission scenarios have shown that medium-to-high navi- 
gation accuracies (40 to 15 nrad in an angular sense) can be achieved through the use of nontraditional 
ata-processing modes for Doppler and ranging data types acquired from the Deep Space Network (DSN) 
[l,2j. These new techniques take advantage of improved calibrations of the limiting ground system error 
sources affecting the data and make use of high-speed workstation computers to reduce the data Studies 
are also being conducted to demonstrate the utility of these alternative data-processing modes with actual 
flight data acquired from the Ulysses and Galileo spacecraft [3,4]. 


These promising research findings prompted this investigation into the use of improved radio metric 
data-processing modes in tracking and navigational support of high Earth orbiter (HEO) missions. In this 
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article, a navigation error covariance analysis is described, which studies the utility of advanced data 
processing modes for X-band (8.4-GHz) Doppler data acquired from DSN-based radio antennas. The 
analysis investigates two-way Doppler-only navigation performance; other radio metric data strategies, 
such as two-way radio ranging, are not addressed. A discussion of Doppler data-processing modes is 
presented along with assumptions for error covariance analysis derived from an improved set of data 
acquisition and orbit determination error-modeling strategies. Results from the covariance analysis are 
described for four sample highly elliptical orbits of the space very long baseline interferometry (SVLBI) 
mission set. A discussion section highlights some critical implications from the analysis and areas that 
will require further investigation. 


II. Doppler Data-Processing Modes 

A. Phase Versus Frequency 

DSN Doppler data acquired in a two-way coherent mode are not direct frequency shift measurements, 
but integral counts of the number of cycles of the transmitted carrier signal relative to the received 
carrier signal that have accumulated over a tracking pass [1]. These cycle counts are differenced to form 
measurements of the average Doppler shift over short time intervals (typically 1 to 10 min); it is these 
differenced- count Doppler measurements that have traditionally been used for navigating spacecraft [5]. 
The purpose of differencing the Doppler counts was to overcome limitations in early tracking hardware 
that caused discontinuities in the data (known as cycle slips ) that occurred when the ground receiver’s 
phase-tracking loop would momentarily lose its lock on the spacecraft carrier signal. Cycle slips are most 
likely to occur when the spacecraft’s Doppler frequency is large and varies rapidly or when the spacecraft 
carrier signal-to-noise ratio approaches the tracking threshold of the ground receiver. 

The notion of using the original Doppler count as a navigation measurement was considered as early 
as 1966 by Curkendall [6], but the idea was not popularized because of the operational constraints 
cited above. Given the steady improvements in the DSN tracking system over the years, together with 
the development of robust methods for resolving occasional cycle slips, there is a renewed interest in 
using the original Doppler count (subsequently referred to as total-count phase) as a viable data type for 
spacecraft navigation [l]. 1 2 The motivation for using total-count phase as a navigation measurement is that 
the precision of these data is very high (a few millimeters at X-band frequencies) whereas differencing 
cycle counts to construct a frequency-formulated observable effectively increases the data noise level. 
Furthermore, unlike Doppler that is differenced-count formulated, phase-formulated Doppler lends itself 
to a simpler data noise model since the errors are uncorrelated. 

B. Filtering Strategies 

The standard orbit determination filtering strategy used by flight project navigation teams treats 
various systematic error sources as unmodeled consider parameters, which are not estimated but whose 
effects are accounted for (i.e., “considered”) in computing the error covariance of filter (estimated) pa- 
rameters [7]. In a consider state analysis, the estimated parameter set s sensitivity to various unmodeled 
consider parameters can be computed via the partial derivatives of the state estimate with respect to the 
consider parameter set [8]. Depending on the magnitude of the resulting sensitivities, the filter-computed 
estimation error covariance is modified to account for unmodeled effects to generate a more realistic es- 
timate of predicted navigation performance. The filter has no knowledge of the unmodeled parameters’ 


1 Historically, the result of differencing cycle counts to form measurements of the average Doppler shift has been referred 
to as differenced -range Doppler, the convention originating from the mathematical formulation of the observable in which 
the station-to-spacecraft range is explicitly differenced over a specified integration, or “count,” time. 

2 T. D. Moyer, “ODE and REGRES Modifications for Processing Block V Receiver Doppler Observables and Total Count 
Phase Observables,” JPL Engineering Memorandum 314-568 (internal document), Jet Propulsion Laboratory, Pasadena, 
California, August 24, 1993. 
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contribution to uncertainty in the state estimate since the modified covariance (the consider covariance), 
including effects from both the estimated and consider parameters, is not fed back to the filter. 

Some principal reasons for using a consider state filter are: (1) certain parameters, such as fiducial 
station locations, may be fixed in order to define a reference frame and/or length scale; (2) there may 
be a lack of adequate models for an actual physical effect; (3) computational limitations exist when at- 
tempting to adjust parameters of high order, such as the coefficients in a gravity field; or (4) if estimated, 
the computed uncertainty in model parameters would be reduced far below the level warranted by model 
accuracy [9,10]. Consider state filters have been known to experience failure modes, such as when addi- 
tional data yield an increase in the consider covariance, or when the consider covariance propagates to an 
unreasonably large result over time [10]. In these instances, it is usually necessary to empirically “tune” 
the filter by adjusting data weights or model assumptions to obtain useful estimates. 

A new sequential data-filtering strategy currently under study is the enhanced orbit determination 
filter, in which most or all of the major systematic ground-system calibration error sources affecting the 
data are treated as filter parameters, along with spacecraft trajectory parameters [1,11,12]. This strategy 
differs from current practice, in which the ground-system calibration error sources are represented as 
unestimated bias or consider parameters and accounted for only when computing the error covariance 
of the filter parameters. The motivation behind the enhanced filter is not so much to improve upon the 
a priori ground system calibrations, but to incorporate a more accurate model of the physical world into 
the filter [1]. 


III. Error Covariance Analysis Assumptions 

Earlier orbit determination studies of DSN-based Doppler tracking for HEO missions focused on using 
the conventional frequency formulation of the Doppler observable along with a standard consider state 
filtering strategy [13]. For this new study, a revised error covariance analysis was performed to quantify 
the navigational utility of a phase formulation of the Doppler observable together with an enhanced data- 
filtering strategy, which treats the principal ground system calibration error sources as filter parameters. 
In this section, assumptions for the revised navigation error covariance analysis are provided, including 
orbit characteristics, data-acquisition and simulation strategies, and filter error modeling. 

A. Sample Orbits 

Four highly elliptical orbits, all derived from the international SVLBI mission set, were selected for 
analysis. The first orbit represents the current operational orbit design for the Japanese MUSES-B 
spacecraft of the VLBI Space Observatory Program (VSOP), targeted to launch in September 1996 
[14]. Nominal orbit parameters are provided in Table 1. The MUSES-B orbit is the lowest of the four 
sample orbits, with perigee and apogee heights of 1,000 and 20,000 km, respectively. The second orbit is 
representative of the latest orbit design for the Russian RadioAstron project [15]. It is the most eccentric 
of the sample orbits, with perigee and apogee heights of 5,000 and 76,800 km, respectively. Table 2 lists 
the nominal orbit parameters for RadioAstron, targeted to launch in 1997 or 1998. 

The two remaining sample orbits are both based on preliminary orbit designs for the Advanced Radio 
Interferometry Between Space and Earth (ARISE) SVLBI mission. ARISE is intended to be a next- 
generation SVLBI mission with a more ambitious set of scientific goals than VSOP and RadioAstron 
[16]. Current ARISE orbit design calls for perigee and apogee heights of about 5,000 and 12,000 to 
50,000 km, respectively; final orbit selection will ultimately depend on the principal scientific objectives 
of the mission. For this analysis, two candidate orbits were assumed with apogee heights of 12,000 
and 50,000 km. Nominal orbit parameters for a mission launch in 2005 are provided in Table 3. Note 
that ARISE will have very stringent orbit determination requirements that exceed the capability of 
current and anticipated ground-based radio tracking strategies, such as two-way Doppler, and will require 
a much more ambitious tracking and navigation strategy. Current design calls for two onboard Global 
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Table 1 . Sample orbit parameters for the Japanese 
VSOP mission (MUSES-B). 


Parameter 

Value 

Nominal launch date 

September 1996 

Initial spacecraft ephemeris 


Semimajor axis 

16,878 km 

Eccentricity 

0.5629 

Inclination 

31.0 deg 

Argument of perigee 

134.24 deg 

Longitude of ascending node 

116.14 deg 

Mean anomaly 

0.0 deg 

Additional parameters 


Perigee height 

1,000 km 

Apogee height 

20,000 km 

Orbit period 

6.06 h 

Table 2. Sample orbit parameters for the Russian 

RadioAstron project. 

Parameter 

Value 

Nominal launch date 

1997/98 

Initial spacecraft ephemeris 


Semimajor axis 

46,778 km 

Eccentricity 

0.7781435 

Inclination 

51.5 deg 

Argument of perigee 

190.0 deg 

Longitude of ascending node 300.0 deg 

Mean anomaly 

0.0 deg 

Additional parameters 


Perigee height 

4,000 km 

Apogee height 

76,800 km 

Orbit period 

28 h 


Positioning System (GPS) receivers [16]. 3 However, for this analysis, it is the orbit characteristics that 
are of principal interest, not the actual mission requirements. 

B. Tracking Data Simulation 

Only two passes of two-way Doppler data were simulated from two different DSN sites. The lengths of 
the data arcs depend on the sample orbit being studied. For shorter orbit periods, such as the MUSES-B 
6.06-h orbit, a 2.25-h pass from the Madrid site and a 2.43-h pass from the Canberra site were assumed. 


3 S. C. Wu and R. P. Malla, “GPS-Based Precision Determination of Highly Elliptical Orbits for Orbiting VLBI Applica- 
tions,” JPL Interoffice Memorandum 335.8-94-004 (internal document), Jet Propulsion Laboratory, Pasadena, California, 
March 29, 1994. 
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Table 3. Sample orbit parameters for the ARISE mission. 


Parameter 

Value 

Nominal launch date 

2005 

Initial spacecraft ephemeris 

Semimajor axis 

18,878.15 km/33,878.15 km 

Eccentricity 

0.40/0.66 

Inclination 

60.0 deg 

Argument of perigee 

0.0 deg 

Longitude of ascending node 

0.0 deg 

Mean anomaly 

0.0 deg 

Additional parameters 

Perigee height 

5,000 km 

Apogee height 

12,000 km/50,000 km 

Orbit period 

7.17 h/17.24 h 


In the case of the ARISE 7.17-h orbit, a 4.07-h pass from Madrid and a 4.20-h pass from Goldstone 
were assumed. The two passes for the 28-h RadioAstron orbit consisted of a 3.03-h Goldstone pass and 
a longer 15.82-h Madrid pass. For the longer 17.24-h ARISE orbit, a 9.52-h Madrid pass and a 4.23-h 
Goldstone pass were assumed. 

To account for random data noise, the measurement error models of [1] were assumed. These models 
are representative of the DSN’s current X-band Doppler system for the 34-m high efficiency (HEF) stations 
with nominal values for spacecraft turnaround ratio, transmit frequency of the carrier signal, and sample 
time. For a differenced-count Doppler measurement at time t k , denoted as f k , the following model is 
assumed: 


fk = pk + Vk 


( 1 ) 


where 


Pk = station-to-spacecraft range rate at time t k 

v k — samples of a zero-mean white Gaussian sequence, in which each sample has constant 
variance and is uncorrelated with all other samples 4 

The random process v incorporates both additive phase- measurement errors and errors due to ground 
system frequency instability that are integrated over the count time of each observation. For this analysis, 
differenced-count Doppler measurements were weighted with a 1 -a measurement uncertainty of 0.1 mm/s 
(metric value) for a 60-s count time. In addition, additive noise variances were adjusted by an elevation- 
dependent function to reduce the weight of the low-elevation data. No data were acquired at elevation 
angles below 10 deg from any DSN site. 

For a total-count phase measurement at time t k , denoted as fa to represent the Doppler count at time 
t k , the following model is assumed: 


This approximation is not rigorously correct since successive differenced-count Doppler data points share common values 
o f the Doppler count and, therefore, each data point is correlated with the two points adjacent to it. In practice, it is 
believed that the uncorrelated measurement error assumption does not yield significantly incorrect statistical calculations 
tor the large Doppler data sets typically used in mission operations [1], 
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0fc = ( Pk ~ Po) + 00 +Vk+€k 


( 2 ) 


where 

p k = station-to-spacecraft range at time tk 
po = station-to-spacecraft range at time to 
eft o = unknown phase offset of Doppler counter at time to 

and 

7 ]k = additive phase- measurement error 
£ k = cumulative phase- measurement error 


The phase offset, 0 O , represents the Doppler counter initialization error and is assumed to be a random 
bias. The r] k samples are assumed to be a white, zero- mean Gaussian sequence with constant variance, 
and the £* values represent the cumulative phase error induced by the integration of frequency variations 
by the Doppler counter. The total-count phase measurements were weighted for this analysis with a 1 -a 
measurement uncertainty of 2.5 mm (metric value). ^ As with differenced-count Doppler, the additive 
noise variances were adjusted by an elevation-dependent function to reduce the weight of low-elevation 
data, and a 10-deg lower elevation cut-off angle assumed for the DSN stations. 


C. Orbit Determination Error Model 

Table 4 provides the dynamic and observational error model assumptions that make up the enhanced 
orbit determination filter, along with a priori statistics, steady-state uncertainties for the Gauss-Markov 
parameters, and noise densities, iV, for the random walk parameters. With the exception of the gravi- 
tational force model, all parameters were treated as filter parameters and grouped into three categories: 
spacecraft epoch state, spacecraft nongravitational force model, and ground system error model. The 
Earth’s gravitational parameter (GM) and geopotential field harmonic coefficients were treated as un- 
modeled consider parameters and grouped in the gravitational force model category. By comparison, 
Table 5 gives the error modeling assumptions that comprise the standard consider state filter model. 


A batch- sequential factorized Kalman filter was used in the estimation process, with a batch size of 
840 min (14 h) for the standard-filtering strategy, reduced to 10-min batch intervals for the enhanced- 
filtering strategy so that short-term fluctuations could be tracked in the transmission media. For process 
noise, first-order Gauss-Markov (exponentially correlated) random processes were assumed. The process 
noise covariance is given by q = (l - m 2 ) a 2 ss where m = exp [- (tj+ i - tj) /r]. Here, tj is the start time 
for the jth batch and r is the associated time constant. The term a S s is the steady-state uncertainty, 
i.e., the noise level that would be reached if the dynamical system were left undisturbed for a time much 
greater than r. For the random walk, both a 93 and r are unbounded (r = oo) and a steady state is never 
reached. The noise density for the random walk is characterized by the rate of change of the process noise 
covariance, q = Aq/At , where A is the batch size and A q is the amount of noise added per batch [9]. 


5 The data noise values for both differenced-count Doppler and total count phase can be readily modified for future navigation 
analyses as the performance specifications of the supporting ground system begin to mature. 
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Table 4. Enhanced orbit determination filter with ground-system error model 
representative of current DSN calibration accuracy. 


Estimated parameter set 

Uncertainty (lcr) 

Remarks 

Spacecraft epoch state 

A priori 

Constant parameters 

Position 

10 3 * * km 

Velocity 

1 km/s 


Nongravitational force model 

Solar radiation pressure 

Steady-state 

Markov parameters 

Specular/diffuse reflectivity 

10 percent of nominal 

0.25-3 day time constant 

Anomalous accelerations 

Steady-state 

Markov parameters 

Radial 

10“ 12 km/s 2 

1-3 day time constant 

Transverse 

10 -12 km/s 2 

1-3 day time constant 

Ground system error model 

Doppler phase offset 

A priori 

Random walk 

(each station) 

100 km 

1 cm 2 /h 

DSN station coordinates 

Crust fixed 

A priori 

Constant parameters 

Spin radius (r fl ) 

0.18 m 


Z-height (z h ) 

0.23 m 


Longitude (A) 

3.6 x 10 -8 rad 


Geocenter offset 

A priori 

Constant parameters 

Z-component 

1 m 


Earth orientation 

Steady-state 

Markov parameters 

Pole orientation 

1.5 x 1CT 8 rad 

1-day time constant 

Rotation period 

0.2 ms 

12-h time constant 

Transmission media 

A priori 

Random walk 

Zenith troposphere 

5 cm 

1 cm 2 / h 


(each station) 


Consider parameter set 

Uncertainty (lcr) 

Remarks 

Gravitational force model 

A priori 

Constant parameters 

Earth’s GM 

GM x 10" 8 


Harmonics 

8x8 field (GEM-L2) 



The principal difference between the enhanced and standard filter models used for this study was the 
modeling of observational errors, namely: 

(1) A random walk model (simple Brownian motion process) was used to track short-term 
fluctuations in the troposphere and assumed zenith delay calibration uncertainties rep- 
resentative of the current DSN-based calibration accuracy. 

(2) A phase offset parameter for each station was included in the ground system error model. 
As with the tropospheric path delays, these parameters were modeled as random walk 
processes. Table 4 shows the noise density given for these parameters, which is derived 
from the white frequency noise representing ground-system frequency instability (see 
Section III.B). 

(3) Three stochastic parameters were included in the ground system error model to account 

for dynamical uncertainties in the Earth’s pole location and rotation period and to rep- 

resent the pole model solutions developed by Finger and Folkner [17]. 
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Table 5. Standard consider-state orbit determination filter with ground-system error 
model representative of current DSN calibration accuracy. 


Estimated parameter set 

Uncertainty (la) 

Remarks 

Spacecraft epoch state 

A priori 

Constant parameters 

Position 

10 3 km 


Velocity 

1 km/s 


Nongravitational force model 

Solar radiation pressure 

Steady-state 

Markov parameters 

Specular/diffuse reflectivity 

10 percent of nominal 

0.25-3 day time constant 

Anomalous accelerations 

Steady-state 

Markov parameters 

Radial 

10“ 12 km/s 2 

1-3 day time constant 

Transverse 

10“ 12 km/s 2 

1-3 day time constant 

Consider parameter set 

Uncertainty (la) 

Remarks 

Ground system error model 

Station coordinates 

A priori 

Constant parameters 

Spin radius (r s ) 

0.18 m 


Z-height ( z h ) 

0.23 m 


Longitude (A) 

3.6 x 10~ 8 rad 


Geocenter offset 

A priori 

Constant parameter 

Z-component 

1 m 


Transmission media 

Zenith troposphere 

A priori 

Constant parameters 

Wet 

4 cm 


Dry 

1 cm 


Gravitational force model 

A priori 

Constant parameters 

Earth’s GM 

GM x 10“ 8 


Harmonics 

8x8 field (GEM-L2) 



(4) The gravitational force model was the same model used in previous studies, with the 
Earth’s GM and truncated (8 x 8) GEM-L2 geopotential field harmonic coefficients 
treated as consider parameters. This is the only element of the overall enhanced fil- 
ter model that used a standard consider state filtering approach. 6 * 


IV. Results 

Results of the numerical error covariance analysis, based on data-acquisition and error-modeling as- 
sumptions described in Section III, are summarized in Table 6. The 1-a position and velocity uncertainties 
for reconstructed orbit estimates are tabulated in a root-sum-square (RSS) sense for two different Doppler 
data-processing modes: (1) differenced-count Doppler data reduced with the standard consider state filter 
and (2) total-count phase data reduced with the enhanced filter. 

The radio navigation performance results in Table 6 assume that science data were collected over 
the radio metric tracking data arcs. Accordingly, it appears that navigation performance is significantly 
improved with the more modern Doppler data-processing mode of total-count phase data reduced with 


6 The argument for not treating the gravitational force model parameters as actual filter parameters is due principally to 

computational limitations when attempting to estimate the harmonic coefficients of the geopotential field (see discussion 

in Section II. B). 
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Table 6. Orbit accuracies of 1-cr for reconstruction over data arc. 


Data-processing mode 

VSOP 

RadioAstron 

ARISE (low) 

ARISE (high) 


RSS position uncertainty, m 


Differenced-count Doppler 
with standard filter 

12.6 

10.9 

11.5 

83.5 

Total-count phase 
with enhanced filter 

5.5 

6.5 

6.5 

30.0 



RSS velocity 

uncertainty, cm/s 

Differenced-count Doppler 
with standard filter 

0.17 

0.03 

0.26 

0.41 

Total-count phase 
with enhanced filter 

0.06 

0.02 

0.13 

0.15 


the enhanced filter. However, these results reflect the accuracies that are achievable only over the specific 
data arcs and not over the entire orbit arcs. More precise representations of the reconstructed orbit 
accuracies over the entire orbit arcs (and in some cases over multiple orbit arcs) are shown in Figs. 1 
through 4. These figures were constructed from filter-generated error covariances, which were smoothed 
and combined with consider parameter sensitivities to produce full consider covariances, then mapped 
orward to give a time history of the reconstructed position and velocity uncertainties over a 24- or 28-h 
period, depending on the sample orbit being evaluated. 

From the time history plots (Figs. 1 through 4), a significant improvement in reconstructed orbit 
accuracies is seen when Doppler data are processed as phase-formulated measurements and reduced with 
t e enhanced filter. This is true for both position and velocity uncertainties for all four sample orbits, with 
the most significant improvement evident in the perigee regions for velocity reconstruction, a regime that 
has historically met with limited success when using traditional radio metric data-processing methods. 

Table 7 attempts to better quantify the performance improvement by giving both the range in uncer- 
tainties and the average (percentage) improvement over the 24- and 28-h time histories. Actual values 
used to generate the percentages of improvement were computed by integrating each error curve over 
the mapped interval to compare total areas of improved versus reference (conventional) data-processing 
modes. From this summary table, relative percentage improvements ranging from about 40 to 60 percent 
are seen, depending on the sample orbit. A slightly more dramatic improvement is seen for reconstructed 
velocity uncertainties over reconstructed positional uncertainties. Again, these results reflect the orbit 
accuracies over the entire propagation or mapping period, i.e., 24 h for the VSOP and both ARISE sample 
orbits and 28 h for the RadioAstron sample orbit. 

Recall that the gravitational force model parameters assumed for both the standard and enhanced 
orbit determination filter models were treated as unmodeled consider parameters (see Tables 4 and 5). 
This was true for both GM and geopotential harmonic coefficients. To gain insight into the effect of 
these consider parameters on the filtering strategy being used, the approximate percentage contribution 
of these error sources on the total reconstructed position and velocity uncertainties over the propagation 
periods was computed, as illustrated in Fig. 5. The altitude dependence of the gravitational force mod- 
eling errors to the total reconstructed orbit accuracies is clearly evident. Not surprisingly, these errors 
contribute far more to the nontraditional data-processing mode; but this is an artifact of the filtering 
strategy being used and not the formulation of the Doppler observable. These results suggest that if im- 
proved navigation accuracies are to be achieved when using the enhanced orbit determination filter, it may 
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Fig. 1. VSOP 1 -a orbit determination accuracy statistics for expected RSS (a) total position 
uncertainty and (b) total velocity uncertainty. 
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Fig. 5. Approximate percentage contribution of 1 -<j gravitational force modeling errors 
(Earth GM and harmonics) to expected RSS (a) total position uncertainty and (b) total 
velocity uncertainty. 




Table 7. Orbit accuracies of 1-cr for reconstruction over entire propagation period. 


Data-processing mode 

VSOP 

RadioAstron 

ARISE (low) 

ARISE (high) 


RSS position uncertainty, m, and 
relative percentage improvement 

Differenced-count Doppler 

4-23 

8-11 

9-60 

24-119 

with standard filter 





Total-count phaser 

2-11, 

4-7, 

4-49, 

8-48, 

with enhanced filter 

56 percent 

43 percent 

41 percent 

63 percent 



RSS velocity uncertainty, cm/s, and 



relative percentage improvement 

Differenced-count Doppler 

0.06-1.8 

0.03-0.5 

0. 2-2.4 

0.3-1. 5 

with standard filter 





Total-count phase 

0.03-0.9, 

0.02-0.2, 

0. 1-2.0, 

0.1-0. 5, 

with enhanced filter 

58 percent 

48 percent 

41 percent 

62 percent 


be necessary to more accurately model the gravitational error sources and possibly treat the relevant 
parameters as actual filter parameters to be estimated along with the spacecraft trajectory parameters. 
The principal motivating factor for using a more sophisticated filtering strategy ultimately depends on 
mission requirements, bearing in mind the altitude dependence of gravitational force modeling errors. 


V. Discussion 

Although the results from two other possible permutations of candidate Doppler data-processing modes 
were not presented— e.g., differenced-count Doppler with enhanced filter and total-count phase with 
standard filter — error covariance calculations performed for these special cases reflect mixed performance 
results. Phase data reduced with the standard filter actually exhibited about a 40-percent worse orbit 
accuracy than traditional frequency-formulated Doppler with the standard filter because the precision 
of these data is very high and, thus, extremely sensitive to unmodeled ground-system calibration errors. 
As studies of interplanetary trajectories have shown, it is necessary to incorporate major ground-system 
calibration errors affecting the data as filter parameters to take full advantage of Doppler phase (without 
artificially deweighting the data). 

When differenced-count Doppler data were used exclusively and reduced with the enhanced filter, there 
was a modest improvement in reconstructed orbit accuracies (~20 percent). However, it is imprudent 
to use this more complicated filtering strategy for very little gain, unless actual orbit determination 
requirements can be easily met. 

Because of the long data arc lengths assumed for the higher orbits considered in this study, concern 
arose as to whether the presence of broken tracking passes might significantly degrade total-count phase 
navigation performance. Therefore, additional error covariance calculations were made for the RadioAs- 
tron orbit case. The phase passes were broken into three shorter intervals of equal length with a 5-min 
break between passes. This resulted in a net loss of about 15 min of data from the original case. An 
independent phase offset parameter was assumed for each pass and the covariance matrix reset at the 
beginning of each track to represent a Doppler count initialization procedure at each station, effectively 
yielding a new phase offset for each pass [1]. Results from this modified tracking scenario exhibited no 
marked degradation in reconstructed orbit accuracies from the original case, despite a 1 5-min reduction 
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in overall data arc length. This result clearly illustrates the robustness of the enhanced filter to solve for 
additional offset parameters incurred from a data-acquisition scenario involving broken tracking passes. 


VI. Conclusions 

A revised navigation error covariance analysis was performed for four highly elliptical Earth orbiters 
derived from the SVLBI mission set. This new study focused on utilizing recently developed or enhanced 
Doppler data-processing modes to reduce X-band Doppler data acquired from DSN-based radio track- 
ing stations. Preliminary error analysis suggests a factor of 2 to 4 improvement in orbit accuracies is 
achievable over traditional data-processing modes when Doppler data are formulated as total-count phase 
measurements rather than differenced-count frequency measurements and processed with an enhanced 
data-filtering strategy that incorporates the major ground-system calibration error sources affecting t e 
data as filter parameters. 

Future work in this area will focus on a thorough sensitivity analysis to determine which dynamic 
and observational sources of error will require further modeling improvement or additional calibration 
accuracy. Plans for concept demonstrations are also being drafted that will use actual DSN-based radio- 
metric tracking data acquired during past mission operations in support of highly elliptical Earth-orbiting 
spacecraft. The JPL operational orbit determination software set is currently undergoing verification and 
validation tests for new upgrades that will facilitate the use of phase-formulated Doppler observables for 
use in both interplanetary and Earth-orbiter mission navigation support. 
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Seamless Data-Rate Change Using Punctured 
Convolutional Codes for Time-Varying 
Signal-to-Noise Ratios 
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In a time-varying signal- to- noise ratio (SNR) environment, symbol rate is often 
changed to maximize data return. However , the symbol-rate change has some un- 
desirable effects, such as changing the transmission bandwidth and perhaps causing 
the receiver symbol loop to lose lock temporarily, thus losing some data. In this 
article, we are proposing an alternate way of varying the data rate without chang- 
ing the symbol rate and, therefore, the transmission bandwidth. The data rate 
change is achieved in a seamless fashion by puncturing the convolutionally encoded 
symbol stream to adapt to the changing SNR environment. We have also derived 
an exact expression to enumerate the number of distinct puncturing patterns. To 
demonstrate this seamless rate-change capability, we searched for good puncturing 
patterns for the Galileo (14,1/4) convolutional code and changed the data rates by 
using the punctured codes to match the Galileo SNR profile of November 9, 1997 . 
We show that this scheme reduces the symbol-rate changes from nine to two and 
provides a comparable data return in a day and a higher symbol SNR during most 
of the day. 


I. Introduction 

In deep-space communications and other space communications, the signal- to- noise ratio (SNR) varies 
during a day. The degree of variation is determined by weather conditions, antenna elevation angle, 
antenna-pointing accuracy (both the transmitter and receiver antennas), changes in satellite latitude, 
and many other factors. For example, the total signal power-to-noise density ratio, P t /N 0 , during a 
typical 24-hour pass for the Galileo Mission can fluctuate in a range between 16 and 22 dB-Hz. In order 
to maximize the data return in this time- varying SNR environment, the transmitted symbol rate is varied 
as a function of the estimated PtjN 0 at the antenna. The symbol rate is set as high as possible under 
the constraint that the symbol SNR is high enough for the tracking loops to remain in lock and that the 
bit-error-rate (BER) requirement is met. In the Galileo Mission, there are six different symbol rates, and 
there can be as many as eight symbol-rate changes (from 10 to 640 symbols/s) during a day. One problem 
associated with the symbol-rate change at a low operating symbol SNR is that the symbol synchronization 
loop may have to reacquire the symbol phase, which may cause real-time data loss. A technique that 
involves opening the symbol loop at the moment of the symbol-rate change has been proposed, 1 but this 
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technique requires very accurate time predicts on the moment of change. It is not clear if the predict 
information can be obtained within the required accuracy. 

In this article, we are proposing a simple and low-cost alternative solution to the data rate-change 
problem by changing the data rate at the error-correcting coding stage rather than at the transmission 
stage. The data rate is changed by puncturing the low-rate convolutional code while the symbol rate is 
kept constant. In this way, the basic structures of the encoder and decoder remain unchanged, making 
the scheme simple and less costly. The idea is to minimize the number of symbol-rate changes and still 
maintain a high enough symbol SNR for the loops to remain in lock and the BER to stay low. Symbol 
rate is changed only if the symbol SNR goes too high (wasting bit SNR) or too low (making the receiver 
unable to track the symbols). 

By allowing the code-rate change, we are essentially adding a degree of freedom in the data return- 
maximization problem. The code rate will share a part of the necessary data-rate changes with the symbol 
rate, therefore reducing the number of symbol-rate changes. This feature becomes even more important 
when the available bandwidth is fixed. 


In Section II, we will present the definition and an overview of puncturing patterns. In Section III, we 
will discuss our procedure for selecting good puncturing patterns, and Section IV will provide an example 
of using the punctured convolutional code for the SNR profile of the Galileo Mission on November 9, 
1997, which is an arbitrarily chosen day. In Section V, we will give some concluding remarks. 


II. Definition and Enumeration of Puncturing Patterns 

A. Definition of Puncturing Patterns 

A regular-rate 1/IV convolutional code generates IV code symbols per bit. By periodically and system- 
atically refraining from transmission of some of the code symbols, a higher rate code can be constructed 
from an original lower rate 1/IV code. Let the period be L bits or NL code symbols. We define a punc- 
turing pattern P of period N L symbols to be an NL binary-tuple, where a 1 denotes that the symbol in 
the corresponding location is to be sent and a 0 denotes that the symbol is to be deleted. If there are m 
zeros in P, the resulting punctured code is a higher rate L/(NL - m) code, where 0 < m < (N - l)L. 
For example, let N = 4, L = 4, and one puncturing pattern be P = {0111 1110 1011 1101}. We define 
the rightmost digit to correspond to the first symbol and the rightmost group of four digits to correspond 
to the four symbols of the first bit. The puncturing pattern, P , indicates that the second symbol in the 
first bit, the third symbol in the second bit, the first symbol in the third bit, and the fourth symbol in the 
fourth bit in a period are not transmitted. The resulting punctured code is a rate 4/(4 x 4 — 4) = 1/3 code. 
With the leftmost digit being the most significant bit and the rightmost digit being the least significant 
bit, the puncturing pattern, P, can be represented as 7 ebd in hexadecimal form. 

B. Enumeration of Puncturing Patterns 

In this section, we develop an exact expression to enumerate the number of unique puncturing patterns. 

Clearly, there are (^ L ) different possible patterns for P. Since P is repeated every L bits or NL 
symbols, any cyclic shift of N symbols in P gives the same code performance as P. However, this does 
/wr Fe( * UCe ^ le num ber of patterns that give a distinct code performance by a factor of IV, as some of the 
( m ) Patterns may have a smaller period L t . That is, L, divides L, which is denoted by L % | L. Let f(L % ) 
denote the number of puncturing patterns with period L l exactly (including 1). Notice that /(L t ) = 0 if 
the m zeros cannot be evenly divided among LjL x partitions (i.e., {L/L t ) /m). Also, among the 
patterns with period L t , some may have smaller periods. Let p be a prime that divides L t . If ^(mP./L), 
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then there are (Jl^Lp)) Patterns of P with period Li/p. The total number of distinct puncturing 
patterns is, therefore, 




Li\L 


where f(Li) can be enumerated as follows: 

HU) = (S) - E (t) 

V L / p \ Li V Lp ' 


Notice that we define the combinatoric function (“) = 0 if either m or n is not an integer. In the above 
example with N = 4, L = 4, and m = 4, an exhaustive search requires checking ( 4 ) = 1820 puncturing 
patterns. By taking into account the cyclic property of the puncturing patterns, the number of distinct 
puncturing patterns is now reduced to 



which is a reduction by almost a factor of 4. 


III. Puncturing Pattern Search Procedure 

In this section, we describe the search procedure that we used to find good puncturing patterns for a 
rate- 1 AN convolutional code. Using this procedure, we searched for punctured patterns for the (14, / ) 
convolutional code used for Galileo. We punctured it from rate 1/4 to rate 1/3, then to rate 1/2. A rate 
compatibility [1] restriction is added in the puncturing-pattern search. That is, a code symbol used in the 
high-rate code is also used in the low-rate code. For example, to search for a rate-1/2 punctured co e, 
we puncture the rate-1/3 code found a step before, not the rate-1/4 code. This was necessary mainly 
because of limited computing resources. 

For each punctured code rate, the goal is to find the puncturing pattern, P, that gives the lowest BER 
at that rate for a range of SNR values. Direct simulation of the punctured convolutional code is not viable 
since there are so many different puncturing patterns. As a first step in selecting the puncturing patterns, 
we computed the weight profile of each punctured code that includes the free distance, df ree , the number 
of paths of weight d, a d , and the information bit error weight, c d . To further simplify our search, we 
only searched for paths of weight d such that d free < d < d cut , where d cut is some predetermined value 
that is large enough to infer the code’s BER performance and small enough to complete the search in a 
reasonable time. Note that there are L different starting points for diverging paths, where L is the period 
of the puncturing patterns. The worst case is considered in comparing the puncturing patterns. 

A systematic search is carried out to find the patterns with the maximum free distance and the 
minimum number of paths of weight d for Viterbi decoding. Three patterns with the largest free distance 
and smallest number of paths of weight d are then simulated to obtain three BER curves. The lowest 
BER curve is selected and used further to compute the BERs for the concatenated codes of convolutional 
code as the inner code and the Reed-Solomon (RS) (255,223) code as the outer code, assuming infinite 
interleaving. 
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Once we have the points of E b /N 0 versus BER of the concatenated code, we fit a curve through these 
points. These curves are used to determine the BER for a given SNR profile of the Galileo Mission. When 
the BER is less than 10 7 , we determine that the code rate can be used in that time period. 

Note that the (14,1/4) Galileo code is used here only to demonstrate the alternative possibility of using 
punctured codes. In fact, the (14,1/4) code is composed of a (11,1/2) convolutional code and the NASA 
standard (7,1/2) code. The NASA standard (7,1/2) code was necessary because the hardware encoder 
on the spacecraft cannot be altered or bypassed. 

A. Upper Bound on Free Distance 

Before searching for the maximum free distances, we compute the upper bounds of the free distances 
to see the effect of the puncturing period on the free-distance bound of the punctured codes. The upper 
bounds on the free distances for convolutional codes can be computed using expressions given in [21. 
Figure 1 shows some of the bounds on the free distances for codes punctured from code (14,1/4), with the 
minimum period from 1 to 4. By minimum period we mean that the period 4 does not include period 1 
or 2. The results show that a shorter puncturing period gives a higher upper bound on the free distance 
but the shorter puncturing period provides a smaller set of possible code rates. 



Fig. 1. Upper bounds on free distance for punctured codes from (14,1/4) code. 

B. Weight Spectra of Punctured Codes 

The parent code in this case has the following polynomials: 2c22, 3d7d, 2 bed, and ldd3, First we 
search for punctured codes from (14,1/4) to (14,1/3) and find the weight spectra corresponding to all 
different punctured patterns. The period in this case is 4, which corresponds to 464 different puncturing 
patterns. We then sort the weight spectra in ascending order according to the number of paths of weight 
, at. Finally, we pick the best three patterns, and their weight spectra are shown in Table 1. According 
^ weight spectra, the best pattern is 6666. This implies that the puncturing pattern has period 1, and 

2cL "d77^d°idd3 UnCtUred ° Ut time ' ThiS COrreSp ° nds t0 the ( 14 ’V3) code with polynomials 
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To further puncture the code to rate 1/2, we use the best 1/3 code found earlier as the parent code. The 
following patterns are found to be the best: 3636, 3535, and 3333 in octal numbers. The weight spectra 
of the three best puncturing patterns are shown in Table 2. Note that when searching for puncturing 
patterns with period 4, those patterns with period 1 and 2 are included. 


Table 1. Weight spectra of punctured codes (14,1/3) from (14,1/4) parent code. 


Pattern 

d 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

bbbb 

a d 

0 

4 

6 

4 

9 

14 

22 

48 

93 

130 

237 

389 

638 

bdbd 

ad 

1 

1 

3 

8 

13 

21 

27 

54 

68 

137 

225 

400 

652 

bbbd 

o-d 

1 

5 

5 

5 

13 

16 

38 

54 

101 

146 

288 

481 

800 

bbbb 

Cd 

0 

14 

18 

18 

55 

72 

122 

322 

641 

920 

1853 

3134 

5530 

bdbd 

C d 

1 

3 

9 

35 

60 

121 

139 

320 

486 

938 

1699 

3150 

5368 

bbbd 

Cd 

5 

14 

19 

24 

77 

91 

240 

347 

724 

1080 

2313 

4067 

7068 


Table 2. Weight spectra of punctured codes (14,1/2) from (14,1/3) parent code. 


Pattern 

d 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

3636 

a d 

0 

2 

6 

10 

24 

51 

142 

344 

824 

1956 

4726 

11363 

3535 

a d 

0 

2 

8 

9 

35 

70 

154 

371 

931 

2286 

5464 

13234 

3333 

a d 

0 

3 

0 

14 

0 

73 

0 

545 

0 

2884 

0 

16679 

3636 

Cd 

0 

5 

20 

70 

146 

354 

1144 

2914 

7780 

20229 

52967 

5525 

3535 

Cd 

0 

9 

37 

53 

251 

550 

1298 

3370 

9353 

25245 

64261 

35749 

3333 

Cd 

0 

9 

0 

71 

0 

520 

0 

4686 

0 

29943 

0 

4011 


C. BER of Punctured Convolutional Code From Simulation 

The weight-spectra search is only the first step in the code puncturing pattern search. To further 
compare their performance, the punctured codes are simulated with an encoder and the Viterbi decoder 
for several bit-SNR values. The traceback length used in the Viterbi decoder in this case is at least 160, 
and the input soft symbols are quantized with 8 bits. The simulated results are shown in Tables 3 and 4. 
Generally, the three puncturing patterns give similar BERs. 

D. BER of Concatenated Code 

Once we obtain the BER from the Viterbi decoder, we can compute the bit-error rate at the output of 
the RS decoder, assuming infinite interleaving using the expression given in [3, p. 256]. In the case of the 
Galileo Mission, there are 8 bits in a codeword, 255 codewords in a frame, and the number of correctable 
errors is 16. The computed BERs at the output of the RS decoder are shown in Tables 5 through 7. 


IV. Example Using the Galileo Profile 

We use the predicted SNR profile of the Galileo Mission on November 9, 1997, as an example to 
explain how the number of symbol-rate changes can be reduced with code-rate changes. For a given SNR 
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Table 3. BER of punctured convolutional codes (14,1/3). 


E h /N„ 

Puncturing patterns 

bbbb 

bdbd 

bbbd 

-1.2494 

0.3528 

0.3532 

0.3533 

-0.7494 

0.2320 

0.2336 

0.2330 

-0.2494 

0.1083 

0.1100 

0.1090 

0.2506 

0.0342 

0.0345 

0.0347 

0.7506 

0.0076 

0.0077 

0.0076 

1.2506 

0.0012 

0.0012 

0.0012 


Table 4. 

BER of punctured convolutional codes (14,1/2). 

E b /N„ 


Puncturing patterns 


3636 

3535 

3333 

-1.0103 

0.4389 

0.4339 

0.4445 

-0.5103 

0.3573 

0.3519 

0.3625 

-0.0103 

0.2230 

0.2206 

0.2279 

0.4897 

0.0924 

0.0934 

0.0924 

0.9897 

0.0230 

0.0243 

0.0226 

1.4897 

0.0040 

0.0043 

0.0035 

1.9897 

0.0005 

0.0005 

0.0004 

2.4897 

4.7 x 10“ 5 

5.4 x 10~ 5 

3.5 x 10 -5 


Table 5. 

BER output of RS decoder using punctured code (14,1/4). 

E b /N a 

BER input to RS decoder 

BER output of RS decoder 

-2.0 

0.4218 

0.4218 

-1.5 

0.3408 

0.3408 

-1.0 

0.2139 

0.2139 

-0.5 

0.1023 

0.1023 

0.0 

0.0326 

0.0194 

0.5 

0.0070 

5.4 x 10" 9 

1.0 

0.0013 

2.3 x 10" 2 ° 

1.5 

0.0002 

1.1 x 10 -49 
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Table 6. BER output of RS decoder using punctured codes (14,1/3). 


E b /N 0 

BER input to RS decoder 

BER output of RS decoder 

-1.2494 

0.3528 

0.3528 

-0.7494 

0.2320 

0.2320 

-0.2494 

0.1083 

0.1083 

0.2506 

0.0342 

0.0230 

0.7506 

0.0076 

1.8 x 10" 8 

1.2506 

0.0012 

1.0 x 10" 20 


Table 7. 

BER output of RS decoder using punctured codes (14,1/2). 

E b /No 

BER input to RS decoder 

BER output of RS decoder 

-1.0103 

0.4389 

0.4389 

-0.5103 

0.3573 

0.3573 

-0.0103 

0.2230 

0.2230 

0.4897 

0.0924 

0.0924 

0.9897 

0.0230 

0.0029 

1.4897 

0.0040 

1.6 x 10 -12 

1.9897 

0.0005 

3.7 x 10 -27 

2.4897 

4.7 X 10' 5 

1.6 x 10“ 44 


profile— for example, the one shown in Fig. 2— the objective is to get the maximum data return under 
the conditions that the bit-error rate is below 1(T 7 and the symbol SNR is maintained above -6 dB for 
the carrier, subcarrier, and symbol loops to track. To achieve this goal, the current plan is to change the 
symbol rate using a fixed code rate, 1/4, and an alternate way is to allow the code rate to change as well, 
thus reducing the number of symbol-rate changes. 

We arbitrarily select a set of three code rates, namely, 1/4, 1/3, and 1/2. The variable code rate can 
only take values from this set. Figure 3 shows the symbol rates using fixed and variable code rates. In 
the fixed-code rate case, there are nine symbol-rate changes, compared to two symbol-rate changes in 
the variable-code rate case. With these symbol rates, each of the two systems will have a symbol SNR 
above -6 dB, as required, where the variable-code rate case has a slightly higher symbol SNR for most 
of the day, as shown in Fig. 4. The code-rate changes are shown in Fig. 5. 

Multiplying the code rates by the symbol rates, we obtain the bit rates as shown in Fig. 6 for the 
fixed- and variable-code rate cases. The areas under the two curves in Fig. 6 are the total data returns 
for the day. The data return using the variable code rate is found to be comparable with that using the 
fixed code rate. 


V. Conclusions 

In this article, we have described a simple and low-cost method to change the data rate to match 
the time- varying P t /N 0 environment. This is done by puncturing the convolutional code at the error- 
correction coding stage rather than by changing the symbol rate at the transmission stage. The main 
advantages of this method are that it allows seamless transition from one data rate to another and that, 
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Fig. 6. Bit rates on November 9, 1997, for Galileo. 


for a fixed available bandwidth, data rate is allowed to change for a larger data return. We applied this 
method to the Galileo SNR profile on November 9, 1997, as an example to demonstrate its effectiveness. 
We showed how this method reduces the number of symbol-rate changes from nine to two and gives a 
comparable data return in a day and a higher symbol SNR for most of the day. 


Notice that, in this example, we arbitrarily picked 100 and 200 symbols/s as two symbol rates to be 
used. We are developing techniques to select the symbol rates that will maximize the data return. The 
problem is formulated below. 


For a given SNR profile, ( P t /N 0 )(t ), we wish to find the symbol rate and the code rate R c {t) 

such that the data return given by R sym (t)R c (t)dt is maximal, subject to the following constraints: 


(1) The number of changes of symbol rate R sym is less than a desired number. 

(2) The BER is below a designed value, BER < BER design . 

(3) The symbol SNR is above the minimum value for the tracking loops to maintain lock, 
E 3 /N 0 > (EJN 0 ) min . 
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Turbo Codes for Deep-Space Communications 

D. Divsalar and F. Pollara 
Communications Systems Research Section 


Turbo codes were recently proposed by Berrou, Glavieux, and Thitimajshima 
[2], and it has been claimed these codes achieve near-Shannon-limit error correction 
performance with relatively simple component codes and large interleavers. A re- 
quired Eb/N 0 of 0.7 dB was reported for a bit error rate of 10~ 5 , using a rate 1/2 
turbo code [2j. However, some important details that are necessary to reproduce 
these results were omitted. This article confirms the accuracy of these claims , and 
presents a complete description of an encoder/decoder pair that could be suitable 
for deep-space applications, where lower rate codes can be used. We describe a new 
simple method for trellis termination, analyze the effect of interleaver choice on the 
weight distribution of the code, and introduce the use of unequal rate component 
codes, which yields better performance. 


I. Introduction 

Turbo codes were recently proposed by Berrou, Glavieux, and Thitimajshima [2] as a remarkable step 
forward in high-gain, low-complexity coding. It has been claimed these codes achieve near-Shannon-limit 
error correction performance with relatively simple component codes and large interleavers. A required 
E b /N 0 of 0.7 dB was reported for a bit error rate (BER) of 1(T 5 , using a rate 1/2 turbo code [2], However, 
some important details that are necessary to reproduce these results were omitted. The purpose of this 
article is to shed some light on the accuracy of these claims and to present a complete description of an 
encoder/decoder pair that could be suitable for deep-space applications, where lower rate codes can be 

US6 | u W ° n<5W contri * 3Utions are reported in this article: a new, simple method for trellis termination 
and the use of unequal component codes, which results in better performance. 


II. Parallel Concatenation of Convolutional Codes 

The codes considered in this article consist of the parallel concatenation of two convolutional codes 
with a random interleaver between the encoders. Figure 1 illustrates a particular example that will 
be used m this article to verily the performance of these codes. The encoder contains two recursive 
inary convolutional encoders, with Mi and M 2 memory cells, respectively. In general, the two compo- 
nent encoders may not be identical. The first component encoder operates directly on the information 
it sequence u = (u i,---,ujv) of length N, producing the two output sequences x lt and xi p . The 
second component encoder operates on a reordered sequence of information bits, u\ produced by an 
interleaver of length N, and outputs the two sequences x 2l and x 2p . The interleaver is a pseudorandom 
block scrambler defined by a permutation of N elements with no repetitions: a complete block is read 
into the ^interleaver and read out in a specified permuted order. Figure 1 shows an example where 
a rate r - 1/n = 1/4 code is generated by two component codes with Mi = M 2 = M = 4, producing the 
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Fig. 1 . Example of an encoder. 


outputs x lt = u, Xip = U • g a /g b , x 2i = u', and x 2p = u' • g a /g b , where the generator polynomials g a and 
g b have an octal representation of 21 and 37, respectively. Note that various code rates can be obtained 
by puncturing the outputs. 


A. Trellis Termination 

We use the encoder in Fig. 1 to generate a ( n{N + M),N) block code. Since the component encoders 
are recursive, it is not sufficient to set the last M information bits to zero in order to drive the encoder 
to the all-zero state, i.e., to terminate the trellis. The termination (tail) sequence depends on the state of 
each component encoder after N bits, which makes it impossible to terminate both component encoders 
with the same M bits. Fortunately, the simple stratagem illustrated in Fig. 2 is sufficient to terminate 
the trellis. Here the switch is in position “A” for the first N clock cycles and is in position “B” for M 
additional cycles, which will flush the encoders with zeros. The decoder does not assume knowledge of 
the M tail bits. The same termination method can be used for unequal rate and memory encoders. 



Fig. 2. Trellis termination. 
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B. Weight Distribution 

In order to estimate the performance of a code, it is necessary to have information about its minimum 
distance, d, weight distribution, or actual code geometry, depending on the accuracy required for the 
bounds or approximations. The example of turbo code shown in Fig. 1 produces two sets of codewords, 
Xl = ( x i», x i P ) and x 2 = (x 2 i, x 2p ), whose weights can be easily computed. The challenge is in finding the 
pairing of codewords from each set, induced by a particular interleaver. Intuitively, we would like to avoid 
pairing low-weight codewords from one encoder with low-weight words from the other encoder. Many such 
pairings can be avoided by proper design of the interleaver. However, if the encoders are not recursive, 
the low-weight codeword generated by the input sequence u = (00 • ■ - 0000100 • • • 000) with a single “1”' 
will always appear again in the second encoder, for any choice of interleaver. This motivates the use of 
recursive encoders, where the key ingredient is the recursiveness and not the fact that the encoders are 
systematic. For our example using a recursive encoder, the input sequence u = (00 • • • 0010000100 ■ ■ • 000) 
generates the minimum weight codeword (weight = 6). If the interleaver does not properly “break” this 
input pattern, the resulting minimum distance will be 12. 


However, the minimum distance is not the most important quantity of the code, except for its asymp- 
totic performance, at very high E b /N 0 . At moderate signal-to-noise ratios (SNRs), the weight distribution 
at the first several possible weights is necessary to compute the code performance. Estimating the com- 
plete weight distribution for a large N is still an open problem for these codes. We have investigated the 
effect of the interleaver on the weight distribution on a small-scale example where N = 16. This yields 
an (80,16) code whose weight distribution can be found by exhaustive enumeration. Some of our results 
are shown in Fig. 3(a), where it is apparent that a good choice of the interleaver can increase the mini- 
mum distance from 12 to 14, and, more importantly, can reduce the count of codewords at low weights. 
Figure 3(a) shows the weight distribution obtained by using no interleaver, a reverse permutation, and 
a 4 x 4 block interleaver, all with d = 12. Better weight distributions fire obtained by the “random” 
permutation {2, 13, 0, 3, 11, 15, 6, 14, 8, 9, 10, 4, 12, 1, 7, 5} with d = 12, and by the best-found permutation 
{12,3, 14, 15, 13, 11, 1,5, 6, 0,9, 7, 4, 2, 10, 8} with d = 14. For comparison, the binomial distribution is also 
shown. The best known (80,16) linear block code has a minimum distance of 28. For an interleaver length 
of N = 1024, we were only able to enumerate all codewords produced by input sequences with weights 
1, 2, and 3. 1 his again confirmed the importance of the interleaver choice for reducing the number of 
low-weight codewords. Better weight distributions were obtained by using “random” permutations than 
by using structured permutations as block or reverse permutations. 



Fig. 3. The (80, 16) code (a) weight distribution and (b) performance. 
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For the (80,16) code using the best-found permutation, we have compared the performance of a 
maximum-likelihood decoder (obtained by simulation) to that of a turbo decoder with 10 iterations, as 
described in Section 3, and to the union bound computed from the weight distribution, as shown in 
Fig. 3(b). As expected, the performance of the turbo decoder is slightly suboptimum. 


III. Turbo Decoding 

Let Uk be a binary random variable taking values in {+1, -1}, representing the sequence of informa- 
tion bits. The maximum a posteriori (MAP) algorithm, summarized in the Appendix, provides the log 
likelihood ratio L(k) given the received symbols y: 


L(k) = log 


P(u fc - +i|y) 
P(u k = — i|y) 


(l) 


The sign of L(k) is an estimate, u k , of u k , and the magnitude |L(fc)| is the reliability of this estimate, as 
suggested in [3]. 

The channel model is shown in Fig. 4, where the nut’s and the n^k's are independent identica lly 
distributed (i.i.d.) zero-mean Gaussian random variables with unit variance, and p = \j2EJN 0 = 
yJ‘2rE h /N n is the SNR. A similar model applies for encoder 2. 

" 1 / 


yi/=py + n i/ 


Vlp = P x 1p + n 1p 


n 1p 

Fig. 4. The channel model. 

Given the turbo code structure in Fig. 1, the optimum decoding rule maximizes either P(ufc|yi,y 2 ) 
(minimum bit-error probability rule) or P(u|yi,y 2 ) (maximum-likelihood sequence rule). Since this 
rule is obviously too complex to compute, we resort to a suboptimum decoding rule [2,3] that sep- 
arately uses the two observations yi and y 2 , as shown in Fig. 5. Each decoder in Fig. 5 computes 
the a posteriori probabilities P{uk\yi, u,), i = 1,2 see Fig. 6(a), or equivalently the log-likelihood ra- 
tio Li(fc) = log(P(u fc = +l|yi,Ui))/(P(«fc = -l|y<,Ui)) where fii is provided by decoder 2 and Q 2 is 
provided by decoder 1 (see Fig. 6(b)). The quantities u, correspond to “new data estimates,” “innova- 
tions,” or “extrinsic information” provided by decoders 1 and 2, which can be used to generate a priori 
probabilities on the information sequence u for branch metric computation in each decoder. 

The question is how to generate the probabilities P(ui,k\uk) that should be used for computa- 
tion of the branch transition probabilities in MAP decoding. It can be shown that the probabilities 
P(ufc|^i,fc) or, equivalently, log (P(uk = +l|iii,fc))/ (P(«fc — — l|«i,fc))> i = 1,2 can be used instead of 
P(«i,fcK) for branch metric computations in the decoders. When decoder 1 generates P(u k \u 2 ,k) or 
log (P(u k - +I|u2,fc)) /(P( u fc = -l|it2,fc)) for decoder 2, this quantity should not include the contribu- 
tion due to which has already been generated by decoder 2. Thus, we should have 



P{u k = +i|u 2 ,fc) _ , P(u/t = + i|yi,«i,i.-- 
og P(u k = -i|« 2 ,fc) og P(«fc = -i|yi.«i,i.- • 


+ ’ * * , Hi,n) 


( 2 ) 
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Fig. 5. The turbo decoder. 
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Fig. 6. Input/output of the MAP decoder: (a) a posteriori probability and (b) log-likeiihood ratio. 

To compute log (P(u k = + I|u 2 ,fc))/ (P(u k = -l|u 2| /c)), we note [see Fig. 6(a)] that 


P(u k |yi,u x ) = 


P\Uk\yuu lyir - • 




^1,1 1 * * * i ui.fc+i, • • • , u \, n ) 


Since u^ k was generated by decoder 2 and deinterleaving is used, this quantity depends only weakly on 
yi an ^ 3 7 ^ h. Thus, we can have the following approximation: 


P(^l,k\ u kyyi,Ux y i, • • • , * 


■ * i ,yv) ~ P(ui jk \u k ) = 2P(u k \u x ^ k )P{u l k ) 


( 4 ) 


Using Eq. (4) in Eq. (3), we obtain 
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P{uk\yi, ui)f’(wi.fc|yi,wi,i> ’ ' ' > «i,fc-i>«i,fc+i> ‘ ‘ ' i“i,w) 
P(u k \y u ui, 1 ,---,u 1<k - U ui,k + i,---,ui,N) - 2P(« fc |ui, fc )i s («i,fc) 

( 5 ) 

It is preferable to work with likelihood ratios to avoid computing probabilities not involving u k (see 
Fig. 6(b)). Define 


Li{k) = log 


P{Uk = +f|M»,fc) 
P(u k - -l|Uj,fc)’ 


i = 1,2 


( 6 ) 


From Eqs. (2) and (5), we obtain V^Hk) = L[ m) (k) - Z/™ 1] {k) at the output of decoder 1, before 
interleaving, for the mth iteration. Similarly, we can obtain L\ m \k) = L 2 \k) — L 2 \k) at the output 
of decoder 2, after deinterleaving. Using the above definitions, the a priori probabilities can be computed 
as 


P(u k = +l|Ui,k) 


1 + e L 'M 


1 - P(u k = -l|ui, fc ), i = l,2 


( 7 ) 


Then the update equation for the mth iteration of the decoder in Fig. 5 becomes 

L[ m) (k) = + a m [j L ( ?\k ) - 4 m) (fc)] . = 1 


( 8 ) 


This looks like the update equation of a steepest descent method, where L^ik) - L[ m \k) 
the rate of change of L(k) for a given u k , and a m is the step size. 


represents 


Figure 7 shows the probability density function of L\(k) at the output of the second decoder in Fig. 1, 
after deinterleaving and given u ^ = H~l. As shown in Fig. 7, this density function shifts to the right as 
the number of iterations, m, increases. The area under each density function to the left of the origin 
represents the BER if decoding stops after m iterations. 

At this point, certain observations can be made. Note that L 2 {k') at the input of decoder 2 includes 
an additive component 2pyhk) which contributes to the branch metric computations in decoder 2 at 
observation y 2 i/c- This improves by 3 dB the SNR of the noisy information symbols at the input of 
decoder 2. Similar arguments hold for L\(k), An apparently more powerful decoding structure can be 
considered, as shown in Fig. 8. However, the performances of the decoding structures in Figs. 8 and 5 are 
equivalent for a large number of iterations (the actual difference is one-half iteration). If the structure in 
Fig. 8 is used, then the log-likelihood ratio L 2 {k) fed to decoder 2 should not depend on Ui* and y[ ik , 
and, similarly, L\{k) should not depend on u 2 k and y 2 i k - Using analogous derivations based on Eqs. (2) 
through (5), we obtain 


L 2 (k) = L 1 (k)-L l (k)-2py' uk 


Li(k) — L 2 (k) — L 2 (k) — 2py' 2ik 


where y' H is the sum of y u with the deinterleaved version of y 2 i and y f 2i is the sum of y 2i with the 
interleaved version of y H . Thus, the net effect of the decoding structure in Fig. 8 is to explicitly pass to 
decoder 2 the information contained in yu (and vice versa), but to remove the identical term from the 
input log- likelihood ratio. 
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IV. Performance 

The performance obtained by turbo decoding the code in Fig. 1 with random permutations of lengths 
N = 4096 and N = 16384 is compared in Fig. 9 to the capacity of a binary-input Gaussian channel for 
rate r = 1/4 and to the performance of a (15,1/4) convolutional code originally developed at JPL for 
the Galileo mission. At BER = 5x 1CT 3 , the turbo code is better than the (15,1/4) code by 0.25 dB for 
N = 4096 and by 0.4 dB for N = 16384. 



So far we have considered only component codes with identical rates, as shown in Fig. 1. Now we 
propose to extend the results to encoders with unequal rates, as shown in Fig. 10. This structure improves 
the performance of the overall rate 1/4 code, as shown in Fig. 9. The gains at BER = 5 x 10 relative to 
the ( 15 , 1 / 4 ) code are 0.55 dB for N = 4096 and 0.7 dB for N = 16384. For both cases, the performance 
is within 1 dB of the Shannon limit at BER = 5 x HT 3 , and the gap narrows to 0.7 dB for N = 16384 
at a low BER. 


V. Conclusions 

We have shown how turbo codes and decoders can be used to improve the coding gain for deep-space 
communications while decreasing the decoding complexity with respect to the large constraint-length 
convolutional codes currently in use. These are just preliminary results that require extensive further 
analysis. In particular, we need to improve our understanding of the influence of the interleaver choice 
on the code performance, to explore the sensitivity of the decoder performance to the precision with 
which we can estimate Eb/N a , and to establish whether there might be a flattening of the performance 
curves at higher E b /N 0) as it appears in one of the curves in Fig. 9. An interesting theoretical question 
is to determine how random these codes can be so as to draw conclusions on their performance based on 
comparison with random coding bounds. 
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In this article, we have explored turbo codes using only two encoders, but similar constructions can 
be used to build multiple-encoder turbo codes and generalize the turbo decoding concept to a truly 
distributed decoding system where each subdecoder works on a piece of the total observation and tentative 
estimates are shared among decoders until an acceptable degree of consensus is reached. 



ENCODER 2 

Fig. 1 0 . The two-rate encoder. 
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Appendix 

The MAP Algorithm 


Let Uk be the information bit associated with the transition from time fc - 1 to time k, and use s as 
an index for the states. The MAP algorithm [1,4] provides the log likelihood given the received symbols 
yk, as shown in Fig. A-l. 


y ^ 

MAP ALGORITHM 

L < u *) 



II 

-A 

£ 


Fig. A-1. The MAP algorithm. 


L(k) = log 


Pjuk = +i|y) 

P{u k = — l|y) 


Ea Ey7+i(yt' s '' s H-i( s ')A( s ) 
g Ea Ea' S'. sK- l(s')0k(s) 


(A-l) 


The estimate of the transmitted bits is then given by si0n[L(fc)] and their reliability by \L{k)\. In order 
to compute Ftp (A-l), we need the forward and backward recursions, 


. . Ea' Ei=±l s j s ) Q! fc-i( s ) 

ak[S) = Ea Ea' 


, , Ea- Ei=±l Ti(yk+l,S, sJ0k+ii/) 
mS) ~ Ea Ea' Ej=±l TjiVk+ly 9 1 , s)ctk(s') 


where 
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7 i{yk,s',s) = ( % e ' 

lo 


■ P J2"=l Vk,uXk,A>' <i) 


if transition s' 
otherwise 


s is allowable for u k = i 


(A-3) 


P = y/2 (EjNo), Vk = P{u k = ±l|u fc ), except for the first iteration in the first decoder, where r/ k = 1/2, 
and x ktl/ are code symbols. The operation of these recursions is shown in Fig. A-2. The evaluation of 
Eq. (A-l) can be organized as follows: 

Step 0: a o (0) = 1 a 0 (s) = 0, Vs / 0 

Pn( 0) = 1 0n(s) = 0,Vs ^ 0 

Step 1: Compute the 7*, ’s using Eq. (A-3) for each received set of symbols y k - 
Step 2: Compute the a fc ’s using Eq. (A-2) for k = 1, • • • , AT. 

Step 3: Use Eq. (A-2) and the results of Steps 1 and 2 to compute the /3*.’s for k - N, ■■ ■ , 1. 
Step 4: Compute L(k) using Eq. (A-l) for k = 1, • • ■ , N. 



Fig. A-2. Forward and backward recursions. 
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for Multiprobe Missions 
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Code division multiple-access spread spectrum has been proposed for use in 
future multiprobe/multispacecraft missions. This article considers a general par- 
allel interference-cancellation scheme that significantly reduces the degradation ef- 
fect of probe (user) interference but with a lesser implementation complexity than 
the maximum-likelihood technique. The scheme operates on the fact that paral- 
lel processing simultaneously removes from each probe (user) the total interference 
produced by the remaining most reliably received probes (users) accessing the chan- 
nel. The parallel processing can be done in multiple stages. The proposed scheme 
uses tentative decision devices with different optimum thresholds at the multiple 
stages to produce the most reliably received data for generation and cancellation of 
probe/spacecraft interference. The one-stage interference cancellation was analyzed 
for two types of tentative decision devices, namely, hard and null zone decisions. 
Simulation results are given for one- and two-stage interference cancellation for 
equal as well as unequal received power probes. 


I. Introduction 

Historically, the Pioneer Venus Mission employed frequency division multiple-access (FDMA) multi- 
plexing to prevent any possible interference between the four probe signals and the two redundant flyby 
spacecraft bus signals. A 2-MHz open-loop recording bandwidth was sufficient to capture all six S-band 
(2.3-GHz) signals. Future missions have been proposed with significantly more probes and use of X-band 
for improved radio metric tracking. The addition of the higher X-band (8.4-GHz) frequency with higher 
Doppler shifts, and the much larger number of probes, 1 implies that a much larger bandwidth would 
be required if an FDMA scheme were used. As a result, Charles D. Edwards proposed that code divi- 
sion multiple-access (CDMA) multiplexing be used for future multiprobe missions. 2 This approach has 
the advantage of greatly reducing the required open-loop recording bandwidth, reducing the cost of the 


1 For example, the Venus Multiprobe Mission (VMPM) is a proposed Discovery Mission concept that delivers 18 small 
probes into the Venus atmosphere and two redundant flyby spacecraft bus signals. The science goal of the mission is to 
understand the superrotation of the Venusian atmosphere, which causes the clouds of Venus to rotate 60 times faster than 
the surface. 

2 C. Edwards, “VMPM Wind Experiment,” viewgraph presentation (internal document), Jet Propulsion Laboratory, 
Pasadena, California, February 24, 1994. 
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recording systems and the complexity of the postencounter data reduction, and simplifying the probe 
design by allowing all probes to use the same ultrastable oscillator frequency and identical transmitter 
structure, with the exception that different seeds should be used for the pseudorandom noise (PN) codes. 

In the following, the term “user” will be used for “probe,” “spacecraft,” or any direct sequence spread- 
spectrum transmitter that communicates with a central receiving Earth station, usually called a “base- 
station receiver.” 

Multiuser communications systems that employ CDMA exhibit a user capacity limit in the sense that 
there exists a maximum number of users that can simultaneously communicate over the channel for a 
specified level of performance per user. This limitation is brought about by the ultimate domination 
of the other user interference over the additive thermal noise. Over the years, researchers have sought 
ways to extend the user capacity of CDMA systems either by employing optimum (maximum-likelihood) 
detection [1) or interference-cancellation methods [2-4], In this article, we discuss a general parallel 
interference-cancellation scheme that significantly reduces the degradation effect of user interference but 
with a lesser implementation complexity than the maximum-likelihood technique. The proposed scheme 
operates on the fact that parallel processing simultaneously removes from each user the total interference 
produced by the remaining reliably received users accessing the channel. In this way, each user in the 
system receives equal treatment in so far as the attempt is made to completely cancel his or her multiple 
user interference. 

When compared with classical CDMA, which has no interference cancellation, and also with the 
successive (serial) interference-cancellation technique previously proposed by Viterbi [3] in which user 
interference is sequentially removed one user at a time (the first user sees all of the interference and the 
last user sees none), the parallel cancellation scheme discussed here achieves a significant improvement 
in performance. Aside from increasing the user capacity, the parallel cancellation scheme has a further 
advantage over the serial cancellation scheme with regard to the delay necessary to fully accomplish 
the interference cancellation for all users in the system. Since in the latter the interference cancellation 

,? nally ’ a deky ° n the ° rder 0f M bit times (M denotes the num ber of simultaneous users in 
the CDMA system) is required, whereas in the former, since the interference cancellation is performed in 
parallel for all users, the delay required is only 1 bit time (for a single-stage scheme). 


II. Single-Stage Interference Cancellation 

A. Tentative Hard Decisions— Equal Power, Synchronous Users 

We consider first the performance of the single-stage parallel interference-cancellation scheme illus- 
trated in Fig. 1, where the tentative decision devices associated with each user are 1-bit quantizers (hard 
ecisions). This particular case corresponds to the scheme proposed in [2] and [4]. We assume that all 
users have the same power; thus, it is sufficient to characterize only the performance of any one user, say 
the first, which will be typical of all the others. Furthermore, we assume that all users have synchronous 
ata streams and purely random PN codes. 3 While the assumption of synchronous users is perhaps un- 
realistic from a practical standpoint, it can be shown that the synchronous user case results in worst-case 
performance and thus serves as a lower bound on the user capacity achievable with this scheme. Alter- 
nately stated, any degree of data asynchronism among the users will yield a better performance, e.g., 
more users capable of being supported for a given amount of signal-to-noise ratio (SNR) degradation’ 
than that arrived at in this section. ’ 

In general, the received signal in Fig. 1 is the sum of M direct-sequence binary phase shift kev (BPSK) 
signals, each with power S u bit time T bl PN chip time T c , and additive white Gaussian noise with 


3 For very long linear feedback shift registers, PN codes can be assumed to be purely random. 
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Fig. 1. A single-stage interference cancellation scheme with parallel processing for CDMA (complex baseband model). 





single-sided power spectral density (PSD) N 0 W/Hz, which at baseband can be written in the complex 


M 

r(t) = £ y/S i m i (t)PN i (t)e’* i + n(t ) 

i=l 


(i) 


where for the ith user, PN t (t) is the PN code; m t (t) = EZ-oo^p(t - kT b ) is the modulation with 
the fcth bit aki taking on equiprobable values ±1 and unit power rectangular pulse shape p(t) of duration 
T b ; and 0, is the carrier phase. For our case of interest here, S, = 5; i = 1, 2, ■ • • , M. After despreading 
and demodulating r(t) with user l’s PN code and carrier reference signal (both of these operations are 
assumed to be ideal), the normalized output of the integrate-and-dump (I&D) circuit is given by 
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where E b — ST b denotes the bit energy; a 0l is the polarity of user V s bit in the interval 0 < t < T b \ n\ = 
^/(y/Tb) / 0 n(t)PNi(t)dt is a zero- mean complex Gaussian random variable with variance E{\n\\ 2 } = No 

representing the thermal noise; and n u = y/S/Tb ft" PN^PN^dt = y/E blu ; i = 2, 3, . . - , M are the 
interference noises contributed by the other M - 1 users, which are modeled as independent zero-mean 
Gaussian random variables, each with variance ST C . 6 Also, the first subscript on x denotes the stage at 
which we are observing the I&D output, while the second subscript denotes the particular user. This 
notation will be useful later on in our discussion of multiple-stage cancellation schemes. The foregoing 
modeling of user interference as additive Gaussian noise follows from the assumptions made in a similar 
analysis of a CDMA system [5], namely, a large spreading ratio, tj = T b /T c , and purely random PN codes. 

Tentative hard decisions are made on the signals, i ~ 1, 2, • • • , Af, and are used in an attempt to 
cancel the other user interference. If a correct tentative decision is made on a particular other user’s 
bit, then the interference from that user can be completely cancelled. On the other hand, if an incorrect 
tentative decision is made, then the interference from that user will be enhanced rather than cancelled. 
A quantitative description of this will be given when we model the signal upon which final decisions are 
made. As we shall see, the performance analysis associated with this model is complicated by the fact 
that the tentative decisions are not independent of one another . More about this shortly. 

After respreading/remodulation, interference cancellation, and despreading/demodulation, the nor- 
malized output of the I&D corresponding to the final decisions is given by 


/i 

S 

M 

Xu = aoi yJ~E b -f- nie ** 1 + ^Eb ^ (3) 

2 = 2 


where 


4 

5 

6 


For convenience, we shall use complex notation to represent the various signals in the receiver. 


Since we are working with a baseband model, the term u remodulation” or “demodulation” refers to complex multiplication 
by the particular user’s carrier phase or its complex conjugate, respectively. 

T -l "°l m /4‘ Zed mterfererlce noises 7i>; » = 2, 3, • • • , M have variance equal to the reciprocal of the spreading ratio, i.e., 
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is a three-valued (0 ±2) indicator random variable whose magnitude represents whether or not a correct 
tentative decision is made on the ith user’s bit. It is tempting to model the ft’s as independent random 
variables Unfortunately, this leads to optimistic results (when compared with the true performance 
results obtained from simulation). In addition to the fact that the ft’s are not themselves independent, 
they are also dependent on the PN cross-correlations, i.e„ the ^s. Fortunately, however, the ft s are not 
strongly dependent, i.e., the only terms that preclude complete independence of, say, ft and ft, are a 0j 
in ft and a 0l 7« = o 0 i7ii in ft. Hence, for sufficiently large M, it is reasonable to assume a Gaussian 
model for the total residual (after cancellation) interference term in Eq. (3). The accuracy of this model 
will improve as M increases (actually, as the number of nonzero terms in /i increases, which implies a 
high tentative-decision error rate). We shall be more detailed about this issue later on when comparing 
the performance results derived from this analytical model with those obtained from a true computer 
simulation of the receiver. 

Assuming then a Gaussian model for h (note that I\ is not zero mean), then the average probability 
of error associated with the final decisions is given by 


P b (E) = -Pr{Re{x n > 0| o 


5 _i}} + ^Pr{Re{x n < O|„ 01 =i}} 


= Pr{Re{x n > 0| aoi= _i}} = Pr|lV t > yf¥ b - y/E b Y, fan cos (ft " <h ) j 
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where 7 


M M 

N t = Re{nie~ j4>l + h - h} = Ni + 0E b y) ft 7ii cos (ft - 4>\) - \/Pb5I^ i7liCOS ^ 1 “ ^ 

1=2 i = 2 

is the effective noise seen by user 1 after cancellation, which in view, of the above, is modeled as a real 
zero-mean Gaussian noise random variable whose thermal noise component AT X has variance o Nl - N 0 / 2. 
It is straightforward to compute the variance of N t as 


= E b (M - 1)0* 7jj cos 2 (0, - 0i) - E b {M - l) 2 ^ft7n cos(0i <f> i)^ 


+ E b {M - 1 )(M - 2)ft7iift7ij cos(ft - 4>\) cos (ft - (pi) (7) 

where i can take on any value from the set 2, 3, ■ • • , M. Hence, from Eq. (5), the average probability can 
be obtained as 

7 To simplify the notation here and in what follows, it is understood that the statistical mean ptfu cos (<j>i - <j>\) is computed 
under the hypothesis aoi = — 1- 


44 


Pb(E) = Q 


(8) 



where 
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“ “ SNR degradation factor (relative to the performance of a single BPSK user transmitting alone) and 
is the Gaussian probability integral defined by 


C(i)4 ^/ exp (- 1 Q <i! ' 


( 10 ) 


Thus, the evaluation of P b (E) reduces to the evaluation of the various statistical averages (moments) 
ot 5 h rec l uired in E q- (9)- These statistical averages, which must be performed over the Gaussian noise 
and interference random variables as well as the uniformly distributed carrier phases, are not trivial to 

C ° m P” te ‘ Nevertheless, they can be obtained in the form of definite integrals of tabulated functions with 
the following results: 
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where ( E b /N 0 ) R denotes the required bit energy-to-noise spectral density ratio for M users communicating 
simultaneously, each of which operates at an average bit error rate P b (E). 

It is common in analyses of CDMA systems [5] to define a degradation factor (loss), D, as the ratio 
(in dB) of the E b /No required to achieve a given bit error rate in the presence of M users, namely, 
( E b /N 0 ) R , to that which would be required to achieve the same level of performance if only a single user 
were communicating, namely, (E b /N 0 )i. By the definition of {E b /N 0 ) u we have 


P b (E) = Q (y/2{E B /N 0 )i) ( 13 ) 

To obtain the degradation factor for a given value of P b (E), we substitute D{E b /N 0 ) i = D 
x \(l/2)[Q- 1 {P b (E))} 2 } for ( E b /N 0 ) R in Eq. (12), which in turn is substituted in Eq. (11). Then, using 
the given value of P b {E), one can solve for D. Unfortunately, a closed-form expression for D cannot be 
obtained, so the results will be obtained numerically. Before presenting these numerical results, however, 
we briefly review the analogous results for conventional CDMA and the successive (serial) interference- 
cancellation scheme proposed by Viterbi [3] (later patented by Dent [6]), since we shall use these as a 
basis of comparison to demonstrate the increased effectiveness of parallel cancellation. 

1. Comparison With Conventional CDMA and Successive Interference Cancellation. In 

a conventional CDMA system, there is no attempt made to cancel the other user interference. Hence, 
(E b /N 0 )i is given by 

/ Eb\ (gbk (E b /N 0 ) R ( 14 ) 

\K) i ~ + (M - l)ST c 1 + (M - l)^(E b /N 0 )R 

Thus, the degradation factor, D, is [5] 


_ ( Eb/Np)^ 1 ( 15 ) 

(E b /No)i 1-{M -Vv-HEb/Noh 


For the successive cancellation scheme [3], Viterbi showed that to guarantee that each user in the 
system sees the same amount of interference from the other users, the user powers should be assigned as 

S k = S l C 1 + ^rT 1 ) , k = M,M~ l,-",2 
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(16) 



where Si is the power of the user to be processed last (the weakest one) and S M is the power of the user 
to be processed first (the strongest one). Distributing the powers as in Eq. (16) ideally guarantees that 
all users see the same ratio of signal power to effective noise spectral density and, thus, the user to be 
processed first (the one that sees all the user interference) is not at any SNR disadvantage relative to the 
user to be processed last (the one for which all interference has been removed). In view of the above the 
degradation factor for the fcth user is given by 


n . (Eb/Noh, _ S* _ /, -l/jr i\t i t * -1 

( E b /N 0 ) x ( E b/ N o)i) (17) 

where ( E b /N 0 ) Rk denotes the required bit energy-to-noise spectral density ratio for the fcth user. The 
average degradation factor, D, for the M user system is obtained by averaging Eq. (17) over k which 
yields 
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M " Mr,~HE b /N 0 )i 


(18) 


It should be emphasized that the result in Eq. (18) ignores the effect of decision errors made at the various 
successive interference-cancellation stages; that is, the interference cancellation is assumed perfect. As a 

result, numerical results derived from Eq. (18) will be optimistic when compared to the actual performance 
of the scheme. 

2. Numerical Results. To illustrate the significant performance advantage of the parallel inter- 
ference-cancellation scheme in Fig. 1 we consider a plot of D versus M for an average bit error probability, 8 
Pb(E) — 10 , and a spreading ratio, r/ = 100. Figure 2 shows the analytical performance of the three 

schemes (conventional, successive interference cancellation, and parallel interference cancellation) as well 
as computer simulation results for the latter. We see that for the conventional and parallel interference- 
cancellation schemes there exists a user capacity limit in that regardless of how much one is willing to 
increase (E b /N 0 ) R (for a given (E b /N 0 ) u or equivalently, a given P b (E)), the required bit error rate 
cannot be achieved if more than M max users simultaneously access the system. For conventional CDMA, 

Mmai ~ 1 + (E b /N 0 )i = 1 + (l/2)[Q-i(P 6 (£))]2 ( 1Q ) 


whereas for the parallel interference-cancellation scheme, the solution is determined from 
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together with the moments in Eq. (11), where now 
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8 The value of Pfc(E) = 10 -2 


is chosen to allow for obtaining computer simulation results in a reasonable amount of time. 
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Fiq. 2. Performance of interference cancellation schemes with equal power users (Viterbi’s scheme with 

unequal power users is also shown). 
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It is emphasized that the user capacity limit for the parallel interference-cancellation scheme comes about 
entirely because of the finite probability of error associated with the tentative decisions. From Fig. 2, it 
appears that the successive interference cancellation does not have a user capacity limit. This is because 
in [3] it was assumed for this scheme that the interference cancellation is perfect, i.e., the effect of decision 
errors at the various interference-cancellation stages was not accounted for. 

Comparing the analytical and simulation results for the parallel interference-cancellation scheme we 
observe that the analytical results are somewhat optimistic. The discrepancy between the two stems 
from the assumption of an analytical Gaussian model for the total residual user interference in Eq. (3), 
whereas the computer simulation makes no such assumption and, thus, predicts the exact performance. ’ 

B. Tentative Hard Decisions— Unequal Power, Synchronous Users 

The results of the previous section can be generalized to the case where the users have unequal powers, 
i.e., Si, i = 1 , 2 , • • • , M. Let ary = Si/Sj denote the ratio of the power of the ith user to that of the jth 
user, who is arbitrarily considered to be the desired user. After interference cancellation, the normalized 
output of the I&D corresponding to the final decisions of user j is, by analogy with Eq. (3), 
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where rij ( \ / f () n(t)PNj(t)dt,j — 1, 2, • • ■ , M is a zero-mean complex Gaussian random variable 
with variance N 0 representing the thermal noise of the jth user; 7 Jt = (1 /T b ) f Th PN J (t.)PN l (t)dt, i ± j 
are the normalized interference noises of the other M - 1 users as seen by user° j (7^ has variance r?" 1 ; 

see Footnote 6), and E bl = S 1 T b is the bit energy of the ith user. Also, analogous to Eq. (4), 3, is now 
defined by 


Pi = aoi - sgn 
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Following steps analogous to Eqs. (5) through (7), we arrive at the desired result for the bit error proba- 
bility of the desired (the j th) user, namely, 


— Q 



(24) 


where (for a 0 j = -1), 
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As an example, consider a group of M users with powers exponentially distributed (linearly distributed 
on a dB scale) over a range of 10 dB between the minimum and the maximum. This model might 
correspond to a distribution of users that are exponentially distant from the base station within a cell. 
Assume that we fix the error probability of the lowest power user (assumed to be user 1 for convenience 
of notation) equal to 1(T 2 (all others would then obviously have a lower error probability). Then, Fig. 3 
illustrates the degradation factor, £>i, of user 1 versus M. For comparison, the results corresponding to 
conventional CDMA with the same user power distribution are also shown in this figure. By comparing 
Fig. 3 with Fig. 2, we observe that in the unequal power case, parallel interference cancellation offers more 
of an advantage over conventional CDMA. The reason behind this observation is that the larger power of 
the other users (which are producing the user interference to user 1) produces tentative decisions with a 
smaller error probability, which in turn results in a better degree of cancellation with regard to the final 
decisions. 


III. Parallel Interference Cancellation Using Null Zone Tentative Decisions 

Much like the idea of including erasures in conventional data detection to eliminate the need for making 
decisions when the SNR is low, one can employ a null zone hard-decision device [see Eq. (27)] for the 
tentative decisions to further improve the fidelity of the interference-cancellation process. The idea here 
is that when a given user’s signal-to-interference ratio is low, it is better not to attempt to cancel the 
interference from that user than to erroneously detect his data bit and, thus, enhance his interference. 
Following the development in Section II. A for a single-stage scheme with equal-power synchronous users, 
then the normalized output of the I&D corresponding to the final decision on user l’s bit a 0 i is still given 
by Eq. (3), with ft now defined by 


ft = a 0 i - nsgn Re 
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\f~Eb a 0i + ^2 

V =5! 



(26) 


where “nsgn” denotes the null zone signum function defined by 


( 1, * > C 

nsgn i = < 0, — C < :r < C 

l-l, x<-< 


(27) 


Here ft takes on possible values (0, ±1, ±2), and its magnitude is an indicator of whether a correct decision 
is made (ith user’s interference is perfectly cancelled), no decision is made (ith user s interference is 
unaltered), or an incorrect decision is made (ith user’s interference is enhanced). Once again, a Gaussian 
assumption is made on the total residual interference; then, since the final decisions are still made as 
hard decisions, the average bit error probability is still given by Eq. (8) together with Eq. (9), with the 
statistical moments of ft, now given by 
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Fig. 3. Performance of interference cancellation schemes with unequal power users. 
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and C = C /VEb is the normalized decision threshold that should be chosen to minimize D for a given 
P b (E) and (E b /N 0 )i determined from Eq. (13). Superimposed on the performance results for the hard 
limiter previously given in Fig. 2 are the results for the null zone limiter. For the specified processing 
gain and average bit error probability, we see that using a null zone limiter allows the maximum number 
of users that can be supported to be increased by about 10 percent. For convenience, the normalized 
threshold has been fixed at C = 0.2. For an unequal (exponentially distributed) power distribution 
among the users, the corresponding results using null zone tentative decisions are superimposed on those 
previously discussed in Fig. 3. For convenience, the normalized threshold has been fixed at C' = 0.4. Here 
again we see a modest improvement in performance. 
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IV. Multiple-Stage Interference Cancellation 


The single-stage scheme of Fig. 1 can be improved upon by cascading multiple stages of parallel 
interference cancellation. The idea here is to repeatedly improve the fidelity of the M tentative decisions 
since each successive stage sees less and less interference. Note that in principle this idea is similar to 
what Viterbi accomplishes in the serial interference-cancellation scheme except that here at each stage 
we simultaneously act on the interference from the most reliable users rather than one user at a time 
An analysis of the performance of such a multistage scheme is difficult if not impossible to obtain due 
o the fact that the tentative decisions at the ith interference-cancellation stage depend on the tentative 
decisions at the ( i - l)st stage. Because of this difficulty, numerical results for the performance of the 
multistage parallel interference scheme will be obtained from computer simulation. Illustrated in Figs 2 
and 3 are performance results for a two-stage parallel interference canceller with hard and null zone 9 
tentative decisions, respectively. We observe that there is significant gain to be achieved by going to more 
than one stage. 


V. Conclusions 

A parallel interference-cancellation scheme was proposed that uses tentative decision devices with dif- 
ferent optimum thresholds at the multiple stages to produce the most reliably received data for generation 
and cancellation of user interference. The one-stage interference cancellation was analyzed for two types of 
entative decision devices, namely, hard and null zone decision. Simulation results are given for one- and 
two-stage interference cancellation for equal as well as unequal power users. The results indicate that, by 
using multiple stages with optimum thresholds at each stage, performance can be significantly improved 
relative to conventional CDMA. Although linear tentative decisions can be used, the performance of such 
a sc h e m e is inferior to one with nonlinear tentative decisions, as our simulations have shown. However 
this scheme with noncoherent detection does not require amplitude and phase estimation. 
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An investigation into the combining of image-processing schemes, specifically an 
image enhancement scheme, with existing compression schemes is discussed. Results 
are presented on the pyramid coding scheme, the subband coding scheme, and 
progressive transmission. Encouraging results are demonstrated for the combination 
of image enhancement and pyramid image coding schemes, especially at low bit 
rates. Adding the enhancement scheme to progressive image transmission allows 
enhanced visual perception at low resolutions. In addition, further processing of 
the transmitted images, such as edge detection schemes, can gain from the added 
image resolution via the enhancement. 


I. Introduction 

There is a new trend developing in the image-processing and image compression fields that has to do 
with the convergence of the two fields. This convergence has now become known as “second generation” 
image coding. It is the result of a growing need to handle large amounts of image data either in transmis- 
sion or in automated image handling— such as image database query and retrieval— while the classical 
compression schemes are reaching their limits. It is now accepted that, in order to achieve more advanced 
compression schemes, we need to use our knowledge about images and their characteristic behavior to 
advantage in compression. 

In this article, we present an initial attempt to combine an image-processing scheme, specifically, an 
image enhancement scheme, with existing image compression schemes. At the image-processing end, 
we use our knowledge about the behavior of edges across scale (across different resolutions) in order to 
extrapolate in scale and increase the resolution of a blurred input image. The ability to extrapolate in 
scale is very useful for compression. We can envision data rate savings by not transmitting at certain 
frequencies and trying to reconstruct the information back at the receiver’s end; we can think about 
combining image processing with existing image compression schemes, such as the subband coding (SBC) 
and pyramid schemes, to achieve additional savings; and, finally, we can use this ability in progressive 
transmission applications, whereby the lower resolution images get enhanced and, thus, information can 
be extracted at earlier stages of the transmission. 

These are some of the issues that we investigate in this article. In Section II, we describe the pyramid 
compression schemes, specifically a variation on the Burt and Adelson scheme [1], and we investigate its 


1 Currently at the California Institute of Technology, Pasadena, California. 

2 Currently with Microsoft Corporation, Redmond, Washington. 
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combination with image enhancement. Section III follows similar steps in relation to the SBC compression 
sc eme. Finally, in Section IV, the combination with progressive transmission for savings in analysis time 
is described. 


II. Combining Image Processing With Pyramid Compression Schemes 

A. The Pyramid Representation 


The pyramid scheme codes an input image in a multiresolution representation via the generation of 
subimages of various scales, as shown in Fig 1 " r " ‘ ’ 

M, 


Here, M Xk [ M Vi denotes subsampling by M Xk and 


yi 


in the x and y directions, respectively. Low-resolution subimages Gk are created by passing G * 
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through a low-pass filter, H , and the decimation box. In the encoder [Fig. 1(a)], we transmit subimages 
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where L k is the difference subimage at the kth level, G k is the low-resolution subimage of the fcth level, 
and G kl is the interpolated version of G k (using filter F). In the decoding part [Fig. 1(b)], we reverse 
Eq. (1) to get the original signal, G 0 . The pyramid representation has been introduced in the literature for 
coding purposes, as it was shown to be a complete representation [1]. Perfect reconstruction is guaranteed 
if there is no quantization of the transmitted data, regardless of the choice of filters H and F. 

B. Compression Via the Pyramid Representation 

The pyramid structure can be used for compression purposes. Using the pyramid coding scheme, 
we decompose the original image into several subimages, with different sizes, and then apply different 
quantization and encoding strategies in the different subimages, depending on the signal characteristics. 
For example, in a linear (e.g., Laplacian) pyramid, the signal variance in different subimages tends to 
be different. Usually, lower-frequency subimages have higher variance. Therefore, we would allocate a 
different number of bits to the subimages (more bits per pixel for higher variance subimages). 

There are several points to be made about compression via the pyramid scheme: 

(1) As stated, this is an oversampled system. This means that the number of output pixels 
at the transmitter is greater than that o f the o riginal image. Specifically, if the original 
size is N x N and the decimation box is |2 ( 2] in every level (K levels overall), then the 
total number of output pixels is 


K 


/ ^ Oi o i 


*=0 


( 2 ) 


Compared to subband and transform coding, which are critically sampled systems, it 
seems that the compression ratios we can obtain using this structure are lower because 
we have to transmit more data. However, due to the following two properties, this is not 
always the case. 
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(b) 



Z-! L 0 


(c) 



(2) No matter how we design the decimation filter H and the interpolation filter F (see 
Fig. 1), we can always obtain perfect reconstruction. Therefore, we can incorporate 
desired filters for H and F such that the signal characteristics in every level are better 
suited for compression. This provides great flexibility. For example, we can apply a very 
long filter or a nonlinear filter for image compression or include a motion-compensated 
filter for video compression. This implies that we can take advantage of some nonlinear 
characteristics of images to help the compression. Compared to subband and transform 
coding, which have a lot of constraints on designing perfect-reconstruction systems, we 
can get many advantages here. 
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(3) In tree-structured subband coding, the quantization noise of the higher level will propa- 
gate to the lower level. This is not desirable, because it is hard to describe the quantiza- 
tion noise behavior if we pass it through too many stages of filters. Especially when the 
quantization stepsize is large or the quantization noise is high, the noise behavior cannot 
be modeled using a simple expression. Therefore, some strange effects become apparent. 
However, in the Laplacian case, we can modify the structure to avoid this problem. As 
shown in Fig. 1(c), we can interpolate the quantized low-passed and decimated signals, 
and then take the difference. The advantage is that we now have exactly the same- 
difference subimages as in the receiving end. Therefore, if we quantize these subimages, 
the quantization noise will not propagate. Overall, it is much easier to control the noise 
behavior. 


C. Gaussian and Laplacian Pyramids 

In this article, we mainly concentrate on Gaussian and Laplacian pyramids, as described by Burt and 
Adelson, commonly referred to as the “Burt pyramid” [1], and Anderson, known as the “filter subtract 
decimate (FSD) pyramid” [2], The Gaussian pyramid consists of low-pass filtered (LPF) versions of the 
input image, with each stage of the pyramid computed by low-pass filtering of the previous stage and 

rDDri P ° nding Subsampl ' ng of the filtered output. The Laplacian pyramid consists of bandpass filtered 
(BPF) versions of the input image, with each stage of the pyramid constructed by the subtraction of two 
corresponding adjacent levels of the Gaussian pyramid. The Burt and Anderson Laplacian pyramids differ 
m the details of when the subsampling step is applied and have slightly different bandpass characteristics. 
I he Burt pyramid follows Fig. 1 exactly, whereas the Anderson pyramid, as defined in Eq. (3), is a 
variation leading to more computational efficiency. In the following, we refer to the input image as 
o, the LPF versions are labeled G\ through G K + i with decreasing resolutions, and the corresponding 
difference images are labeled L 0 through L K , respectively. A recursive procedure allows for the creation 
of the Anderson pyramid, as follows: 


G° n+1 = W*G n 

L n = G n — G° +1 (3) 

Gn+i = Subsampled G° + 1 

where G n is termed the nth-level Gaussian image and L n is termed the nth-level Laplacian image. Gen- 
erally, the weighting function, W , is Gaussian in shape and normalized to have the sum of its coefficients 
equal to 1. The values used for the LPF, which is a 5-sample separable filter, are (1/16, 1/4, 3/8, 1/4, 
1/16). Figure 2 presents an example of a Laplacian pyramid representation. 



Fi p. 2 - «dge maps Presented from left to right are the Laplacian pyramid components 

Lq, L, , and t- 2 , respectively. The pyramid components have been appropriately expanded to match in size. 
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It has been shown [1] that the Laplacian pyramid forms an overcomplete representation of the image, 
thus enabling full reconstruction. The reconstruction process entails adding to a given LPF version of 
the image, Gjv, the bandpass images, L n (n = N - 1, ■ • • ,0), thus reconstructing the Gaussian pyramid, 
level by level, up to the original input image, G 0 . This is a recursive process, as in Eq. (4): 

G n = L n + G( n +i)n n — N — 1, * * • , 0 (4) 


where G( n +i)j is the interpolated version of G n +i* 

D. How Does Image Enhancement Come In? 

Several things can be noted from the above description of the pyramid representation. First, we note 
that the Laplacian pyramid consists of the edge maps of the input image at the different resolutions (see 
Fig. 2). We also note that when coding the image into its Laplacian components, most bits need to be 
allocated to the L 0 level. The first observation leads us to the idea of combining our knowledge about 
edge behavior across scale. The second observation allows the combination with compression. Our goal 
is to learn about the behavior of edges across scale so that we can “predict” the L 0 level of the pyramid 
using only the lower-resolution edge maps up to L\. 

1. The Image Enhancement Scheme. In our image enhancement work [3], we concentrate on the 
edge representation of an image across different image resolutions. Edges are an important characteristic 
of images, since they correspond to object boundaries or to changes in surface orientation or material 
properties. An edge can be characterized by a local peak in the first derivative of the image brightness 
function or by a zero in the second derivative, the so called zero crossings (ZCs). An ideal edge (a step 
function) is scale invariant in that no matter how much one increases the resolution, the edge appears the 
same (i.e., remains a step function). This property provides a means for identifying edges and a method 
for enhancing real edges. 

We concentrate on the edge representation of an image across different image resolutions. For this we 
view the image in a multiresolution framework via the Gaussian and Laplacian pyramids. 1 he Laplacian 
pyramid preserves the shape and phase of the edge maps across scale (see Fig. 2). 

The application of the Laplacian transform to an ideal edge transition results in a series of self-similar 
transient structures, as illustrated in Fig. 3(a). An edge of finite resolution would produce a decrease in 
amplitude of these transients with increasing spatial frequency, with the magnitude of the edge going to 
zero at frequencies above the Nyquist limit [see Fig. 3(b)]. An edge of finite resolution can be created by 
starting with a low-resolution Gaussian image and then adding on all the bandpass transient structures. 
To create an edge with twice the resolution requires the creation of a self-similar transient at the next 
level, hereby referred to as j. The most essential features of these transient structures are that they 
are of the same sign at the same position in space; hence, their ZCs line up, and they all have roughly 
the same amplitude. The precise shape of the structures need not necessarily be maintained so long as 
their scaled spatial frequency responses are similar. The simple procedure described next creates localized 
transients for L_i that satisfy all these constraints except for the maintenance of constant amplitude. 
While more complicated procedures could handle the amplitude constraint, it was found that sharpening 
the stronger value edges produces in itself visually pleasing results. 

The pyramid representation can be viewed as a discrete version of the scale-space description of ZC 
that has been introduced in the literature [4-6]. The scale-space formalism gives the position of the ZC 
across a continuum of scales. One of the main theorems states that ZCs of an image filtered through 
a Gaussian filter have nice scaling properties, one of which is that ZCs are not created as the scale in- 
creases. If an edge appears at lower resolutions of the image, it will consistently appear as we shift to 
higher resolutions (see Fig. 3). Although theoretically defined, not much work has yet taken advantage 
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(a) (b) 



Fig. 3. Laplacian transform on (a) an ideal edge and (b) an edge of finite resolution. 


of the image representation across scale. In our work, we utilize the shape invariant properties of edges 
across scale based on the pyramid representation and in agreement with the consistency characteristic of 
the scale-space formalism. 

The objective is to form the next higher harmonic of the given signal while maintaining phase. Figure 4 
illustrates a one-dimensional high-contrast edge scenario. The given input, Go, is shown in (0) of the 
figure, together with its pyramid components, L 0 and G u shown in (1) and (2), respectively. From 
the pyramid reconstruction process, we know that adding the high-frequency component Lq to the G\ 
component can sharpen G\ to produce the input Go- Ideally, we would like to take this a step further. We 
would like to predict a higher-frequency component, L_j, preserving the shape and phase of L 0i as shown 
in (3) of the figure, so that we can use the reconstruction process to produce an even sharper edge, which 
is closer to the ideal-edge objective, as shown in (4) of Fig. 4. The L_\ component cannot be created by 
a linear operation on the given L 0 component (i.e., the frequency spectrum cannot be augmented using a 
linear operator). We can, thus, never hope to create a higher-frequency output by a linear enhancement 
technique. 

It remains to show how the L-\ component of the pyramid can be generated. We extrapolate to 
the new resolution by preserving the Laplacian- filtering waveform shape, together with sharpening via 
a nonlinear operator. The waveform as in (5) of Fig. 4 is the result of clipping the L 0 component, 
multiplying the resultant waveform by a constant, q, and then removing the low frequencies present (via 
bandpass filtering) in order to extract a high-frequency response. 

Equation (5) formalizes the generation of L_ 


L-i =a{C{L 0 )) (5) 

where G(S) is defined as 


( T if S' > T 
C{S) = l S if -T < S < T 
l -T if S < — T 

Here, T = OM(G 0 ) max . 
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(2) PYRAMIDAL COMPONENT, G-| 


Fig. 4. The one-dimensional ideal-edge scenario. 

Generating the new output image entails taking the L_! image as the high-frequency component of the 
pyramid representation. Based on the reconstruction capability of the pyramid representation [Eq. (4)j, 
the new output is generated next as the sum of the given input, G 0 and L_i, as in Eq. (6): 


Enhanced Image = G - 1 = L-i + Go 


(6) 


2. Enhancement Results. We next show experimental results that indicate that the enhancement 
scheme augments the frequency content of an input image, achieving a visually enhanced output. 

A rock scene example is displayed in Fig. 5. Figure 5(a) presents the enhancement results; Fig. 5(b) 
displays the corresponding power spectral characteristics. The blurred input, which can be the result of 
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Fifl* 5. A rock scene example: (a) enhancement results and 
(t>) corresponding power-spectrum characteristics. In each of 
the above figures, the blurred input and original image are 
presented (top left and top right, respectively), followed by the 
enhanced output (bottom). Both visual perception enhancement 
and power-spectrum augmentation are evident. 


cutting off high frequencies due to bandwidth considerations or a “zoom-in” application, is presented at 
e op- eft corner The original image, which we are assuming is not available to the system and which 
we wish to reproduce, is presented at the top right. The result of applying the previously presented 
algorithm to the blurred input is depicted in the bottom of each figure. We get an overall enhancement 
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perception. The enhanced image very closely matches the original one, and the power spectrum of the 
enhanced image is very close to the original power spectrum. For additional enhancement results, the 
reader is referred to [3,10]. 

In conclusion, the enhancement scheme addresses the most important features (edges) required in 
producing enhanced-resolution versions of existing images. The simplicity of the computations involved 
and ease of implementation enable it to be incorporated into real-time applications. 

E. Combining Image Enhancement With Pyramid Coding 

In this section, we combine the image enhancement scheme, described above, with the pyramid coding 
scheme. We have shown the possibility of predicting the L 0 level of the Laplacian pyramid using lower- 
resolution edge maps. The next step is to code an image with and without the L 0 component and evaluate 
the corresponding rate-distortion performance, i.e., investigate the compression savings versus the output 
image results that we can achieve. 

We decompose the original image, Go, into {L 0 , L\, L 2 , G 3 }. We scalar quantize L 2 , L\, and L 0 
and then compute the entropy of the quantized signals. For G 3 , we first apply differential pulse-code 
modulation (DPCM), then compute the entropy. The average entropy of these subimages represents the 
rate (bits per pixel). Here we use peak signal-to-noise ratio (PSNR) as the distortion criterion, defined 

as 


PSNR = 10 log 10 — — 3 P-P ( 7 ) 

Ar j-ij-i 

where I l3 is the ijth pixel of image /, and X and Y represent the horizontal and vertical size of the input 
image, respectively. 

The Lenna image is used for this coding task. Figure 6 presents the rate-distortion curves for Lenna. 
We concentrate first on the two pyramid-coding curves. We note that in using the L 0 component we 
have all pyramid levels and, thus, the reconstruction would be exact apart from the quantization errors 
induced. Using a predicted L 0 (i.e., the actual L 0 component is not being used), we introduce additional 
noise in the reconstruction process. In general, we note the slow degradation of the rate-distortion curve 
using the enhancement scheme, as opposed to the almost linear drop of the original (nonenhanced) curve. 
Of even more interest is that, at very low bit rates, the ability to estimate the L 0 component from the 
given L x component, or the ability to extrapolate in frequency space, allows for better PSNR. 

An example of two images, with and without the Lq component (top left and top right, respectively), 
is shown in Fig. 7, as compared with the original Lenna image (bottom). Both images are coded with 
approximately 1 bit/pixel. We have 0.99 bit/pixel with PSNR = 34.77 dB for the enhanced image and 
1.053 bits/pixel with PSNR — 33.54 dB for the image decompressed with all its L t components. In this 
case, we get a better PSNR and better perceived similarity to the original for the enhanced image with 
the predicted L 0 component than for the image with all components present. This is a very interesting 
and encouraging result. 

Next, we compare the pyramid compression with the discrete cosine transform (DCT) (refer again to 
Fig. 6). The DCT clearly “wins” the PSNR comparison. We note that, at the very low bit rates, the 
differences are quite minimal. In addition, we need to compare the actual images, as opposed to the PSNR 
ratios, as is shown in Fig. 8. We note that the blockiness with the DCT is very evident and possibly more 
distracting to the eye than the artifacts introduced by the pyramid-plus-enhancement scheme. A zoom-in 
image taken from Fig. 8 is presented in Fig. 9. In the DCT coding scheme, we can see strong blocking 
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Fig. 6. Rate distortion curves for the Lenna image. 

effects in the quantized image (this is the case especially when the bit rate is low or when we zoom the 
image up). This phenomenon results from the independent quantization of blocks. This reinforces the 
claim that the PSNR does not in all instances match our visual perception. 

A similar investigation is done on a moon image, whose rate distortion curves are shown in Fig. 10. 
As before, we note the slow degradation of the rate-distortion curve with the enhancement. We see again 
that, at low bit rates, better PSNR is achieved by predicting the L$ component via the enhancement 
processing stage. When comparing it with the DCT rate-distortion curve [Fig. 10(b)], we notice that at 
very low bit rates we actually achieve better performance than the DCT. 

We conclude this section with a few of the moon images. Figure 11 displays the slow degradation 
phenomenon. Two images are displayed: The left one has 1.27 bits/pixel with PSNR = 31.69 dB, 
and the right image has almost half the bit rate, at 0.65 bit/pixel, with a very similar PSNR value of 
31.1 dB. The two images look identical. In Fig. 12, we compare the pyramid scheme (left) to DCT 
(right) at the low bit rate of 0.47 bit/pixel. The PSNR ratio is larger for the pyramid coding in this 
case, PSNR = 30.53 dB, whereas the DCT case has PSNR = 30.49 dB. The blockiness of the DCT 
is certainly visible here. Note the blocks on the main rocks, which actually degrade the possibility of 
identifying rock boundaries, etc. 

In this section, we have shown encouraging results in combining the image enhancement scheme with 
pyramid coding. The results are interesting especially at low bit rates. Possible modifications to the 
pyramid structure can help us in achieving higher compression ratios. Different filters can be applied in 
the pyramid structure to get better performance. We know that no matter what we put in the decimation 
and interpolation filters (see Fig. 1), perfect reconstruction is always guaranteed if there is no quantization 
and transmission loss. This provides great flexibility, since we can design a better filter to achieve better 
rate-distortion performance however the distortion is defined. For example, we can apply a nonlinear 
filter to take advantage of the nonlinear features of the human visual system. We can then obtain better 
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Fig. 7. Comparison of the Lenna image with and without Lq. The top left image includes Lq in the compression, 
the top right uses a predicted Lq, and the bottom is the original Lenna image. 


image Quality with a greater perceptual effect. The use of a median filter has been shown to achieve such 
an improvement [7], The combination of the enhancement scheme with modified pyramid coding schemes 
remains to be investigated. 


III. SBC Schemes 

A. Introduction 

SBC schemes have recently aroused much attention in the areas of image and video compression. There 
are several advantages to these coding schemes that make the technique attractive. Recently, the Motion 
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Fig. 12. Comparison of pyramid (left) versus DCT (right) at low bit rates. 


Pictures Experts Group (MPEG) has proposed using this technique for its audio compression part. It is 
the belief of many people that SBC may take over DCT as a new image and video compression standard. 
Generally speaking, the compression capability of SBC is fairly good. Among all the linear multiscale 
techniques— such as DCT, the Laplacian pyramid, and SBC— it has been shown that SBC can provide the 
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best compression ratio in the rate-distortion mean-squared-error (MSE) sense, although the computational 
complexity is typically higher than that of DCT in achieving this. The picture quality generated by this 
technique is good compared to the annoying blocky effect generated by DCT. One whole frame is processed 
and quantized at a time as opposed to the block-by-block DCT process. However, there is another kind o 
distortion based on aliasing effects that degrades the picture quality substantially when the compression 
ratio is high or the bit rate is low. This is due to the signal loss in the subbands, which results in the 
aliasing cancellation effect provided by the filter bank being lost. This phenomenon gets more visible 
when the quantization noise is higher. It is generally agreed that the picture quality provided by SB 
is better than that of block coding techniques. An additional advantage of SBC is that it is intrinsically 
progressive. This is achieved by dividing the original images into subimages in different frequency bands. 
One then sequentially transmits subimages with, usually, increasing frequencies. Progressive transmission 
is a desired property for many applications, such as data browsing and image frame conversion between 
different signal formats like high-definition television (HDTV) and standard TV. Although DCT can be 
implemented progressively, it is not as immediate a process. 


The basic principle of SBC is, like other linear techniques, to take advantage of the nonuniform 
distribution of the signals’ power spectrum. It is well known that the power spectrum of image signals 
tends to be nonuniform. We first use a filter bank (and decimators) to decompose the original image into 
several subimages in different frequency bands. One then allocates different numbers of bits in different 
bands depending on the signal variance in that band, thus achieving compression. It can be shown that 
in the uniform filter bank case, the quantization step size in every band has to be equal in order to 
obtain minimum MSE. However, taking account of the features of the human visual system (HVS), the 
quantization step sizes in different frequency bands should be different. It has been shown in testing that 
HVS is more sensitive to the lower frequency components. Therefore, we usually set the step sizes in 
lower bands to be smaller. Although not always true, typically, longer tapped filters can provide better 
rate-distortion performance due to better energy compaction. 


There is much literature available on the principles, implementation techniques, and different types of 
multirate filter banks and SBC [8]. We use octave wavelet (tree-structured wavelet)-based decomposition. 
As shown in Fig. 13, we decompose the low-low subimage into four subimages of the next level, i.e., 


LLi — > LLx+i , HLi +\ , LHi+i , HH{+\ 


by a specific filter (e.g., for the HL subimage, it is a high-pass filter in the x direction and a low-pass filter 
in the y direction) and a decimator that decimates the signals by 2 in both x and y directions. Therefore, 
if we decompose the image up to N levels (level 0, ■ • • , N - 1), then there will be 3V + 1 subimages. At 
the receiving end, we recursively reconstruct the signals by 


LLi = £(LL X + 1) + £(LH i+1 ) + £(HL X + i) 4- £(HH l+1 ) 


where £() represents expanding the signals first (inserting a zero between two pixels), then passing them 
through the corresponding filter that was originally used in the decomposition. 
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Fig. 13. SBC block diagram. 

B. The Combination With Image Processing 

As for the pyramid concepts of the previous section, we wish to explore the possibility of using signal 
correlations across frequency bands to be able to extrapolate from a lower frequency band to a higher 
frequency one, thus eliminating the need of transmitting all the bands. 

T* 16 main difference between this case and the Laplacian pyramid case is that here we have the three 
1 erence subimages {LH,, HL U HH t ) per scale instead of the single bandpass (BP) L % component. We 
propose to estimate a higher-frequency level from a lower-frequency one by expanding the low-frequency 
evel (using £) and adding the different subimage components, as in the SBC reconstruction scheme, with 
an additional intermediate enhancement step. 

Several possibilities come to mind in this new scenario: 

(1) We can estimate each higher-level, i, subimage component from its corresponding com- 
ponent in the previous level, i -\- 1, as follows: 


LH l+l — ► LH U HL l+l — ► HLi, HH l+1 — > HH , 


where LH X represents the low-high (x -y directions) difference subimage in level i. 

(2) We can combine all BP components in level (i + 1) and then estimate the corresponding 
subimage summation in level (z): 


LHi+i + HLi+i + HH t +i — ► LH l + HL t + HH t 


(3) We can use only the low-low components and try to predict their behavior across the 
levels, i.e., LL l+l — ► LL { . 
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Continuing work on the Lenna image, Fig. 14 (top row) displays level 2 of Lenna in the SBC de- 
composition. Figure 14 (center row) shows the estimation of the level-1 subimages, following expansion 
and enhancement of the corresponding level-2 images, while the bottom row presents the original level-1 
subimages. There is both a resemblance and a difference between the estimated bandpass images and the 
original ones We wish to see how well we can estimate level-0 images based on the estimated level-1 com- 
ponents, i.e., using information from components in level 2 alone. Adding level-2 components together 
we generate an exact LL X image (see Section III. A). We next expand that image and each of the estimated 
level-1 BP subimages of Fig. 14 (center row) and sum these together to generate the estimated LL 0 image 
[option (1) above] as displayed in Fig. 15 (left). Aliasing effects are evident. Another experiment is to 
take the expanded LL\ and enhance it [option (3)]. This produces the LL 0 image at the right of Fig 15. 
Again very strong aliasing effects are evident. In order to first add all the bandpass components together 
and then expand [as in option (2)], we need to first define an appropriate filter for the combined subbands. 
The £ function of the original SBC decomposition uses specific filters for the specific subbands and, thus, 

is not suitable for the task. 

In the remaining experiments, we examine the possibility of expanding a given LL X image via Gaussian 
interpolation followed by the enhancement procedure. Figure 16 presents the results of such a procedure 
on Lenna. The expanded image is presented in the top left, and the enhanced image is shown m the top 
right. This result very closely matches the original LL 0 image, which is displayed at the bottom of the 
figure We have thus achieved a good reconstruction of the 0-level image based on the information in leve 
2 alone, eliminating the need to transmit the level-1 subimages. Quite surprisingly, the similarity we see 
between the generated and original images is not reflected in the PSNR ratios. Considering the Lenna 
image the PSNR value for the blurred input (top left) as compared to the bottom image is 2o.l dB. 
PSNR for the enhanced image (top right) is 24.58 dB. This is quite unexpected, especially as we look at 
the global statistics of the three images as presented in Table 1. We note that the enhanced image as 
statistical characteristics that closely match the original input image. 

We can conclude the following: 


(1) The PSNR estimate relates to local pixel-value discrepancies between two images. Since 
the enhancement scheme does not attempt to reconstruct back exact pixel values of the 
original image, we are well aware that it is not ideal for the PSNR measure. Still, it is 
very apparent that the PSNR also does not represent the human perception (as in the 
example of Fig. 16). It is quite clear with this example that the estimated (enhanced) 
image is very close to the original. Because of this mismatch, we choose not to produce 
rate-distortion curves with the PSNR measure. 

(2) Overall we can conclude that there is no immediate step to be taken from the pyramid 
representation case (Section II) to the SBC case. This is probably most evident due to 
the aliasing effects. Ignoring certain bandpass components seems to be more crucial since 
these components are needed for dealiasing. All SBC subimages are needed in order to 
eliminate the aliasing (a theoretical analysis of this characteristic can be found in [9j). 
We can compensate for the missing frequencies, but unless more work is done on how to 
resolve the aliasing problem, this issue remains the main obstacle. 

The overall conclusion at this time is that the particular image-processing scheme suggested in this 
article is not applicable, without major modifications, to the SBC scheme. 


IV. Image Enhancement and Progressive Transmission 

We conclude this article with an additional aspect of possible interest, which is the combination of 
image- processing with progressive image transmission schemes. Here, instead of looking for additional 
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Fig. 14. SBC decomposition of Lenna. Difference subimages of levels 2, estimated level 1 
and original level 1, top to bottom, respectively. 


Table 1. Global statistics of the images in Fig. 16. 


Characteristic 

Original 

Expanded image 

Expanded plus 
enhanced image 

Mean 

1.344 x 10 2 

1.344 x 10 2 

1.344 x 10 2 

Standard 

deviation 

41.3 

39 

40 

Power 

140.6 

140 

140.3 
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Fig. 15. Prediction level 0 of Lenna — aliasing effects. 

bit compression capabilities, we are interested in achieving successive approximations in time. By this 
we are referring to progressively transmitting information, from low resolution to high resolution, with 
the desire to extract information during the transmission without waiting to receive the high-resolution 
image. Moreover, we would like to determine at an early stage of the transmission process if the image 
is at* all of interest, so as to determine if the high-resolution image is to be transmitted. 

In Fig. 17, we demonstrate the combination of the integer subband coding (ISBC) scheme 3 of the 
Gaspra image with the enhancement scheme. We note the possibility of detecting the craters and other 
points of interest much more clearly in the enhanced images, even at extreme compression ratios. For 
the scientist, this can be a tool for determining interest in the region. If it does look interesting, the 
full-resolution image can be transmitted, without any loss. 

Another domain of interest is detection of objects in a given scenery. For this task, an initial phase of 
edge detection is usually performed. We have looked into the combination of an edge detection scheme 
with the enhancement scheme to allow for better and quicker object detection. 

A. Combining Edge Detection and Image Enhancement 

The purpose of combining edge detection and image enhancement is two-fold: First, this scheme can 
be applied to progressive multiresolution image-compression systems (and progressive transmission), such 
as the (integer) subband coding schemes and the Laplacian-based pyramid coding scheme, to detect the 
edges of low-resolution images received at the early stages of transmission. We need an edge detector to 
catch the locations of the desired objects as soon as the low-resolution images are available. Scientists 
can then select and send back only the images (and acquired resolutions) of interest, thus saving in the 
required transmission power. 

The second purpose is that, after combining the edge detection and image enhancement schemes, we 
get more enhanced and clear images as compared with the original images. By doing so, we can capture 
many more details in the received images and get a highly detailed edge map. 

3 K.-M. Cheung, “Low-Complexity Progressive Image Transmission Schemes for Space Applications, .JPL Interoffice Mem- 

o rand um 331-93.2-064 (internal document). Jet Propulsion Laboratory, Pasadena, California. 
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Fig. 16. Predicting level 0 of Lenna image from level 2 information: expanding LL 1 level (top left) 
with enhancement (top right) as compared to the original image (bottom). 

An edge is usually defined as the point where a significant change (normally, intensity) occurs In 
the results presented, we use a gradient-based method. Here, several difference (high-passed) filters are 
es,g ned at 4 5 -deg orientation interval preferences. After passing the original images through these filters, 
we add all the filtered (absolute) values in every direction and take a threshold. If the added sum is larger 
than the threshold, then the corresponding pixel is claimed as an edge point. This scheme achieves results 
comparable to common edge-detection schemes found in the literature. Additional details can be found 


We apply the edge detection scheme on a low-resolution image: Figure 18(a) shows a low-resolution 
image of an oilfield image, and Fig. 18(b) shows the detected edge map. It is expected that fewer details 
(e ges) can be captured, due to lost high-frequency components. However, we can improve this by first 
passing the low-resolution image through the image enhancer and then performing the edge detection, 
r igure 18(c) shows the enhanced image, and Fig. 18(d) shows the detected edge map. We see that more 
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Fig. 17. Combining image enhancement with progressive transmission. 


Enh 
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Fig. 18. Edge detection of the oilfield image: (a) a low-resolution image; (b) the detected edge map for (a); 

(c) an enhanced image; and (d) the detected edge map for (c). 

details can indeed be captured. For example, we can count more oil tanks (top), detect the ships more 
clearly (center), and perceive more of the building structures (bottom). 

V. Summary and Conclusions 

In this article, we have done a preliminary analysis on the combination of image compression with 
image-processing schemes, specifically image enhancement. Encouraging results have been achieved with 
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the pyramid compression scheme, especially at low bit rates (which is the new frontier in the compression 
field). The combination of the image enhancement scheme and SBC was not as successful. Overall, 
we conclude that a scheme that needs all frequency bands for dealiasing cannot be easily altered and, 
specifically, cannot have bands removed and processed without introducing aliasing effects. Finally, the 
case for including enhancement in progressive image transmission was made, with results indicating the 
enhanced visual perception at low resolutions. In addition, further processing of the transmitted images, 
such as basic edge detection schemes, can gain from the added image resolution via the enhancement. 

Several issues for further exploration stem from this work. In the pyramid scheme, we are interested 
in the possibility of pursuing steps similar to the ones described, at lower resolutions, i.e., predicting one 
level from a lower- resolution one, at the low resolutions of the pyramid, thus extending the compression 
from the L 0 level to L\, etc. Initial investigation indicates that this is not a simple extension to the 
existing algorithm. As the resolution is decreased substantially, it is more difficult to locate the edges. In 
addition, the sharpening process will require more investigation as to how to “fill-in the regions in the 
image that have been blurred and now have been sharpened. Overall, this idea requires further research. 

From the SBC investigation, we learn that it could be the case that the compression and image- 
processing schemes cannot be combined as is. It is not always possible to take the existing algorithms 
and combine; rather, we might need to rethink the compression scheme together with the image-processing 
algorithms to generate new compression schemes. 

The work presented is very preliminary work. Still, we believe that the results are interesting enough 
to support future work in this direction. We have only touched upon one category of image-processing 
schemes, image enhancement. Other processing, such as actual segmentation of images based on content, 
“model-based coding,” and more, is attracting much interest in the research community as the new 
frontier for image compression. 
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Reflector antenna optimization schemes using array feeds have been used to re- 
cover antenna losses resulting from antenna distortions and aberrations and to gen- 
erate contour coverage patterns. Historically, these optimizations have been carried 
out using the antenna far- held scattered patterns. The far-held patterns must be 
calculated separately for each of the array feed elements. For large or complex 
antennas (which include beam-waveguide antennas), the far-held calculation times 
can be prohibitive. This article presents a method with which the optimization can 
be carried out in the antenna focal region, where the scattering calculation needs 
only to be done once independently of the number of the elements in the array. 
This article also includes the results of a study, utilizing this unique technique, to 
deter mi ne the capabilities and limitations of using array feeds to compensate for 
gravitational induced losses in large rehector antennas. 


I. Introduction 

Array feeds for reflectors have a number of important uses, which include (1) generating contour 
coverage patterns, (2) correcting for reflector distortions, and (3) improving wide-angle scan. Typical 
methods for optimizing the array feed, for each of these applications, are very efficient when a fixed- 
array geometry is utilized and only the feed excitation coefficients are optimized. For this case, only one 
calculation of the radiating fields for each array element is required. For example, to maximize gain in a 
given direction, the optimization can be as simple as taking the complex conjugate of the secondary fields 
resulting from the illumination of the reflector in the given direction by each of the array feed elements. 
For most existing methods, an optimization that allowed the element spacing and size to vary would be 
extremely time consuming since a radiation integral evaluation would be required for each feed element 
at each step of the optimization process. 

A new method of computing the performance of reflector antennas with array feeds is presented that 
obviates the need to recompute the reflector radiation fields when the feed element size or spacing is 
varied. This allows the optimization technique to efficiently include size and spacing as parameters. 

The mathematical formulation is based upon the use of the Lorentz reciprocity theorem, which con- 
volves the focal-plane field distribution of the reflector system with the feed-element aperture field dis- 
tribution to obtain the element response. Thus, the time-consuming reflector-system radiation integral 
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evaluation is only done once for a given scan direction or reflector surface distortion for all array feed 
geometries and types considered. Table 1 tabulates the relative merits of focal region analysis (sometimes 
referred to as reverse or receive scattering) and far-field analysis (sometimes referred to as forward or 
transmit scattering) A rough rule of thumb is that any analysis that requires evaluating the radiated 
fields at a number of far-field observation points is best done with conventional far-held analysis How- 
ever, analysis that requires investigating many different feed or array feed designs is accomplished faster 
and more easily by using focal plane analysis. 

Examples are given using the technique to design an array feed for the correction of gravity-induced 
istortions of a large dual-shaped ground antenna, both conventional and beam waveguide (BWG) and 
to design an array feed for improved wide-angle scan. The main emphasis of the examples, however, is on 
beam-waveguide designs as exemplified by the Deep Space Network (DSN) Deep Space Station (DSS)-13 
34-m antenna at Ka-band (33.67 GHz). 


Table 1. Comparison of antenna computational techniques 
for a single reflector design. 


Characteristic 

Forward method, 
far field 

New reversed method, 
focal plane 

Fields calculated 

Far field 

Focal plane 

Number of far- 
field points 

Any number 

One 

Number of focal- 
field points 

Not applicable 

Any number 

Number of individual 
feeds 

One 

Any number 

Number of array 
combinations 

One 

Any number 

Typical computation 
time for 34-m BWG 
antenna at 33.67 GHz 
(Cray Y-MP2), a h 

6-12 

6-12 

Time to calculate focal 
plane currents, h 

Not applicable 

3-4 

Time for additional far- 
field points 

Trivial 

6-12 

Time for additional feeds 
or array cases, h 

6-12 

0. 1-0.5 


a Due to priority restrictions for large jobs, turnaround times can be 
more significant than actual computational times. 


II. Focal Plane Analysis 

The calculation of the gain of an antenna system, using the antenna focal plane analysis technique 
can be illustrated by a simple single- reflector antenna, as shown in Fig. 1. Referring to Fig. 1(a) the 
first phase of the analysis consists of several steps. First, the focal-plane fields produced by a plane wave 
impinging upon the reflector antenna aperture are computed. Second, the aperture fields of the feed horns 
located at the focal plane are determined, and these fields are then convolved with the focal plane fields 
to provide a set of complex feed weights. The process can be explained as follows: Consider a reflector 
antenna fed by a horn. We wish to determine the gain of this system in a given direction, (6 0 , <p« ) in the 
receive mode. First, consider the Lorentz reciprocity theorem: 
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In this expression, E a and H a are fields radiated by a set of sources J a and M a , and also E h and H b are 
fields radiated by sources J b and M b . The left integral is over a closed surface that encloses the volume 
defined by the integral on the right side. Over an infinite region, the surface integral becomes zero. 1 he 
Lorentz reciprocity theorem can, therefore, be rewritten as 


JJJ {E a J b -H a ■ M b )dv = JJj {E b • J a H b ■ M a) dv (2) 



which says that, as long as the relationship between the fields and their source currents holds true the 
results will_be the same wherever the integrals are evaluated in the region. Let us redefine the a sources 
as sources J ha and Mha, to be associated with the feed horn apertures and generating the fields ~E ha and 
H ha in the aperture plane of the feed. In turn, let us redefine the b sources as sources ~J lp and M< to 
be associated with the antenna reflector system when illuminated by an incident plane wave source from 
a direction (<9 0 , <t> 0 ) and evaluated in the reflector system focal plane. The feed horn apertures are defined 
to be coplanar with the antenna focal plane. 

Since the integration is limited to the aperture plane/focal plane, the integrals reduce to surface 
integrals. Each integral is proportional to the feed-horn output voltage [1], The left-hand equation is 
used in this work since the program that generates the focal-plane equivalent currents outputs currents 
and the program that computes the feed-horn aperture distributions outputs fields. The expression 
relating the feed-horn outputs, v TA , to the currents from the antenna reflector system and the feed- horn 
aperture fields is then 


Vr A oc {E ha ■ j fp - Hha ■ M j v }ds (3) 


The E ha and Hj, a should_be determined in the presence of the antenna reflector system, and the focal 
plane currents Jf p and M/ p should be obtained when the feed horn is present. Such computations 
would require that the interactions between the feed horns and the reflectors be taken into consideration. 
Taking into account these interactions seriously complicates the analysis and increases the computational 
time. Often E ha and H ha are approximated to the aperture fields^ of a horn radiating into an infinite 
homogeneous free space (no reflector), and the focal-plane currents J fp and M fp of the antenna reflector 
system are also obtained in the absence of a feed. This is a reasonable assumption when the feed and 
antenna reflector system are widely separated in terms of wavelengths. 

As will be seen later, there is also a need to obtain the performance of a feed horn in the presence of 
a plane-wave incident field arriving from a direction (6 p ,<p p ) and in the absence of the antenna reflector 

system, Fig. 1(b). In the same manner as presented above, it can be shown that the output voltage for 
such a feed horn is 


v 


rp 



H ha ' 


(4) 


where J pw and M pw are the currents in the feed aperture plane due to the incident plane wave field. 

In Eqs. (3) and (4), the proportionality constants should be the same, being a function of the horn 
aperture characteristics. The proportionality constants can be eliminated by performing the ratio of 
Eqs. (3) and (4) as follows: 


v 


r A 


V 


rp 


ff{Eha * J fp Hha ‘ M f p }ds 

s 

If {Hha ' J pw — Hha ' M p W }ds 

3 


(5) 


Let us now consider two transmit situations. First, let the feed horn radiate in the absence of the antenna 
reflectors and let the radiated field at (r,O p ,<j> p ) be E h and the power be P Q . Then the gain of the feed 
horn in the direction (0 pj 0 p ) is 
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(6) 


Gh 


4irr 2 

r]P 0 


\Eh? 


Next, let the horn illuminate the reflector system. The scattered field at (r, 6 0 ,<p 0 ) is E a . Assume that 
the power that is radiated by the feed horn is still P Q . Then the gain of the complete antenna system in 
the direction (do, 0 o) is 


G 


a 


4irr 2 

rjPo 


\Ea \ 2 


and, consequently, 


G a = G h 



(7) 


( 8 ) 


Prom reciprocity, we know that the radiated fields and horn output voltages are related by 

1 u t~a | 2 / Q \ 

\E h \ 2 " \v Tp \ 2 

Therefore, by combining Eqs. (5) and (9), the overall gain of the reflector antenna system can be found 
in the receive mode from 


G a =G h - 


| ff {Eha ' J fp — Eha ' M f p }ds\ 


| ff {E ha ■ Jpw ~ H ha ■ M pw }ds\ 


( 10 ) 


It should be noted that (0 o ,0o) and (d p ,<p p ) need not be the same. Therefore, (9 p ,<p p ) has been set to 
(0.0, 0.0) for simplicity of analysis when evaluating the interactions of the array feed with an incident 
plane wave and computing the feed horn far-field gain, Gh- 


III. Optimization Technique 

The optimization technique used in obtaining the maximum gain for an antenna reflector system 
illuminated by a group or an array of feed elements is referred to as the conjugate weight or match method 
[2]. Consider a reflector antenna with an array feed operating in the transmit mode and with each feed 
element excited with equal signal levels. In some direction (0 O ,0 o), in which the maximum output fields 
are desired, the output field for each feed element is determined. Let /, represent the complex output 
field voltage of the antenna for feed element i in the direction (0 o ,<t>o)- Then the maximum or optimum 
output field in that direction would be 


v t = £/;/» ( u ) 

i=l 

If in the receive mode, c, is the output voltage from the ith feed element of the reflector antenna system, 
when illuminated with a plane wave arriving from the direction (9o,<Po)', then the total received signal 
from the antenna where each feed element is weighted by its complex conjugate is 
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( 12 ) 


N 

V r = ][> t 
1=1 

From reciprocity, Cij f % = const for all i. Therefore, except for a constant, the expressions that are based 
on the complex weighting in the receive mode are identical to those in the transmit mode. Thus, v T 
also represents an optimum gain solution. In this analysis, the effects of mutual coupling between array 
elements have been ignored. For the size and type of feed elements considered in this study, this is not a 
limitation. 

If the integral portion of Eqs. (3) and (4) are rewritten as follows, 


{E h a t ■ J f p - H hat ■ M f p }ds (13) 


{ ^ ha, ' J pw E ha t ’ M p W }(ts (l^) 


by the use of Eq. (12), Eq. (10) can be rewritten as 


N 

E 

i~ 1 

I £ c *i d i I 

i = l 

This expression is used to determine the optimum antenna gain simply by knowing the focal plane currents 
of the antenna, the array feed geometry, and the feed element aperture fields. The gain of the array feed 
in the absence of the antenna reflectors, G is obtained by first performing a physical optics integration 
over the E ha fields in each feed element aperture to obtain the far-field pattern for each element. Then 
the power and peak fields are computed from the total fields from all elements in the conventional manner 
and used to compute Gh ■ 

The analysis method consisted of computing the focal plane currents of a reflector system using reverse 
(receive-mode) scattering programs. Physical optics is used for a single-reflector antenna design and a 
combination of geometrical optics off of the main reflector and physical optics off of a subreflector is 
used in a dual-reflector antenna system. Both electric and magnetic focal plane currents are computed 
on a fixed grid over which the array-feed element aperture distributions are superimposed. This grid 
becomes the integration grid over which the convolution of the focal plane currents and feed aperture 
field distributions are integrated. The feed aperture fields are in effect interpolated to fall on the points 
established by the reverse scattering program’s focal plane grid. To generate the feed element aperture 
fields at the required grid locations, the far field element patterns are expanded into a set of circular 
waveguide modes. These modes can then be evaluated at the required grid locations to obtain the feed 
aperture fields. 


IV. Accuracy Considerations 

To determine the accuracy of the focal plane technique, the gain of a single reflector antenna was 
calculated using three different approaches: (1) the focal plane technique combined with a physical optics 
program that took a plane wave impinging on the reflector and computed the required focal plane currents, 
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(2) a forward scattering program, using a physical optics technique that made use of the Jacobi/Bessel 
series to describe the output fields [3], and (3) a forward scattering program, using a physical optics 
method that used a triangular integration grid [4]. The gain pattern was calculated for a number o ar- 
field observation points. Obviously, for the forward approach, this simply meant specifying the locations 
of the observation points. For the focal plane approach, however, calculating the gain pattern required 
the incoming plane wave direction to be adjusted such that the direction of propagation was from t e 
observation point. Table 2 summarizes the results. As can be seen, all three techniques agreed within 
a few hundredths of a dB for pattern gain variations of 7 dB. The only precautions needed to get this 
accuracy were to make sure sufficient integration points were used and that all series representations of 

fields had converged. 

Table 2. Accuracy of the focal plane method, antenna gain. 


Pattern 
angle, deg 

Classical forward method 
Jacobi/Bessel, dB Standard PO, dB 

Reversed method 
focal plane, dB 

0.00 

44.16 

44.15 

44.17 

0.35 

43.15 

43.14 

43.16 

0.60 

41.14 

41.14 

41.15 

0.90 

37.14 

37.14 

37.15 


V. Use of Optimization to Minimize Beam Scan Loss 

One application of the optimization approach described in this article is to minimize the scan losses of 
an antenna where the antenna beam has been scanned off the axis of symmetry by a lateral displacemen 
of the feed system. Depending on the size of the displacement, the losses can be quite large, due to the 
optical aberrations produced. The test case used consisted of a 423.55 cm-diameter parabolic reflec or 
with a focal length-to-diameter ratio (F/D) of 0.5. The compensating array feed consisted of up to 
37 elements. The main beam was scanned 6.15 deg off the axis of symmetry to produce a scan loss on the 
order of 8 dB The primary array element that produces the specified beam scan was located 24.36 cm 
off the antenna axis. The antenna parameters and the improvement in performance due to array feed 
optimization is shown in Table 3. The gain for several numbers of array elements is shown. The gain 
for an antenna with the single-feed element is 44.8 dB. With complex conjugated weights applied to 
37 elements, the gain is increased to 51.9 dB for a net improvement of 7.1 dB. 


VI. Compensating Reflector Distortions by Optimization 

As was mentioned in the introduction, the main purpose of the focal plane optimization technique 
was to reduce the amount of computer time required to optimize the design of an antenna system to 
maximize the gain. This is particularly true for studying large, complex antennas such as the 34-m beam- 
waveguide dual-shaped reflector antennas, operating at 33.67 GHz, that are used at the Jet Propulsion 
Laboratory/NASA deep space tracking network. However, a second advantage of the focal analysis is 
that it allows the actual focal field distributions to be displayed. This can many times give some insight 
into why the antenna behaves the way it does. In the conventional far-field optimization approach [5], 
a full scattering calculation is required for each array feed element, for each antenna configuration, and 
for each array geometry. In the focal plane analysis, a scattering calculation is required only for each 
antenna configuration. Any number of array geometries and any number of feed elements can be studied 
without any further scattering calculations. 

Although considerable effort was expended to develop the focal plane analysis technique, the main 
purpose of the study was to determine the properties of a shaped antenna illuminated by an array feed 
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Table 3. Compensating for beam scan loss using 
conjugate weights. 


Number of elements 
in array 

Pattern 
gain, dB 

Gain 

improvement, dB 

1 

44.8 

0.0 

7 

49.1 

4.3 

19 

50.4 

5.6 

37 

51.9 

7.1 

Geometry 

Reflector diameter 

423.55 cm 


F/D 0.5 

Array element diameter 1.59 cm 


Frequency 11.8 GHz 


Array center offset 

24.36 cm 


Beam peak direction 

6.15 deg 



where the array feed’s purpose is to compensate for various antenna aberrations or distortions. The focal 
plane technique is simply the enabling technology. A conjugate weighting technique is used to determine 
the excitation weights for each array feed element so that the antenna losses might be recovered. The 
goal was to compensate the antenna for gain losses that resulted from gravity-induced distortions as 
a function of the antenna’s elevation angle, where the surface is adjusted for minimum distortion at a 
45-deg elevation angle. The main application is for a beam-waveguide (BWG) antenna where the array 
feed system is mounted in a room below the antenna. However, to expedite the validation of the focal 
plane analysis technique, the initial calculations were limited to the main focal plane of the antenna, i e 
no BWG. Figure 2 illustrates the geometry of the DSS-13 34-m BWG antenna, the properties of which 
were made the basis of this study. The frequency used was 33.67 GHz. Of particular interest is the focal 
region at FI, where the initial calculations were made. A second set of calculations was made at F3, the 
ocal region created by the beam- waveguide system in the basement or pedestal room of the antenna. 
Although the beam-waveguide system consists of two parabolic mirrors, one elliptical mirror, and three 
plane mirrors, only the three curved mirrors are included in the analysis. The three plane mirrors, if 
properly sized, should not affect the performance of the antenna system. Figure 3 shows an array feed 
with the RF front-end setup at F3 that was used in another study for an experimental evaluation of 
using array feeds to compensate for antenna distortion losses. Figure 4 shows a seven-element array with 
4.45-cm diameter dual-mode horns, used in the experimental system. This illustrates a typical application 
tor the results of the trade-off study discussed later in this article. 

A summary of results at Fl and F3 for 7.5- and 45-deg elevation angles is shown in Table 4. The 
results are based on the use of optimum feed horn sizes and focus. Both the focal plane or reversed 
scattering method and the forward or conventional scattering method were used. The results for the 
forward method are from an optimization study done at the time DSS 13 was completed. As can be seen, 
the results for the forward and reversed approaches are very close. 

A. Optimization at FI 

The case used for the evaluation consisted of an array feed of seven equally sized elements in a circular 
cluster on a triangular grid. The first step was to select an element size that minimized the antenna loss 
at one of the distortion extremes, such as at an elevation angle of 7.5 deg. A set of 13 element diameters 
was selected, ranging from 2 to 13 cm. For each element diameter, the focal plane optimization technique 
was used to determine the feed element weights that gave the best performance improvement. Figure 5 
contains two curves. The first curve shows the performance of the antenna that results from a single 
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Fig. 3. An array feed with the RF front-end set up at F3. 



GAIN, dB 






Table 4. Gain comparisons for a shaped antenna between BWG and 
non-BWG designs using focal plane or far-field analysis. 21 


Location 

Calculation 

method 

Evaluation 
angle, deg 

One 

element, dB 

Seven 

elements, dB 

Antenna focal point FI 

Reverse 

45.0 

81.120 

81.129 

Antenna focal point FI 

Forward 

45.0 

81.109 

— 

BWG focal point F3 

Reverse 

45.0 

81.039 

81.046 

BWG focal point F3 

Forward 

45.0 

81.044 

— 

Antenna focal point FI 

Reverse 

7.50 

79.442 

79.733 

BWG focal point F3 

Reverse 

7.50 

79.338 

79.557 

a Antenna surface shape i 

is optimized at a 45.0-deg elevation angle. 


on-axis element as a function of element size. The second curve shows the best performance for an array 
of equally sized elements. From the curve, it was found that an array of elements with diameters of 
9.26 cm gave the best performance at an elevation angle of 7.5 deg. Calculations using the optimization 
technique were repeated for a series of antenna elevation angles, using the 9.26 cm-diameter elements, 
to determine the best performance that can be obtained with the optimum element size determined at 
7.5 deg. Figure 6 summarizes the result. The first curve illustrates the performance that would be 
expected due to reflector distortions as a function of elevation angles if a single feed horn were used. The 
second curve shows the improvement in performance that can be obtained with the seven-element array 
using the optimum element weights computed by the optimization technique. The maximum improvement 
at an elevation angle of 7.5 deg is 0.29 dB out of a distortion loss of 1.68 dB. It was found that if the 
number of elements in the ring around the center element was increased from 6 to about 12 elements, an 
additional improvement of 0.18 dB could be obtained for an overall improvement of 0.47 dB (Fig. 7). 

As indicated earlier, the availability of the focal plane fields from the focal plane analysis program 
can give an insight into why the antenna performs the way it does. Figure 8 shows the focal plane field 
distribution at FI for an elevation angle of 45 deg. Overlaying the field distribution are the outlines of 
the optimum-size feed horns. As can be seen, the center horn covers the field distribution from in excess 
of 5.0 dB down to almost -15.0 dB. The outer ring of feed horns covers a region of ripples indicated by 
the multiple contour rings at -22.5 and -30.0 dB. These multiple rings also indicate that phase reversals 
are present in the region of the outer horns. Thus, the center horn captures the majority of the focal 
plane fields, whereas the outer horns couple very little of the fields because of the low field strengths and 
the fact that the phase reversals will cancel out a large part of what field levels are available. Hence, the 
outer elements will not improve the antenna performance, as can be seen in Fig. 6. Figure 9 is a similar 
plot for an elevation angle of 7.5 deg. In this case, the center horn captures less energy since the peak 
level has dropped and the contour associated with the edge of the horn aperture is now only -5.0 dB. 
Although the outer ring of horns covers a region of levels primarily between -5.0 and -25.0 dB, which 
appears to represent the majority of the energy missed by the center horn, as is seen in Fig. 6, only a 
small part of the lost energy was recovered at 7.5 deg. This is due to the fact that the typical feed-horn 
aperture distribution is Gaussian and couples the most energy when associated with a Gaussian focal 
plane distribution, such as is seen in Fig. 8 for the center horn for the 45-deg elevation case. In the case 
of the outer horns at 7.5 deg, the focal plane distribution they see is essentially wedge shaped in a radial 
direction and, therefore, these horns are very inefficient in coupling the energy within their apertures. 

A question raised during the study was whether the array feed could recover losses due to diffraction 
from the subreflector support tripod. Calculations were made to address this question at an eleva- 
tion angle of 45 deg. To simplify the calculation and to avoid developing a special program to compute 
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Fig. 8. Focal plane field distribution overlaid with the 
outline of the 9.258-cm diameter horns at FI for a 45-deg 
elevation angle. 
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Fig. 9. Focal plane field distribution overlaid with the 
outline of the 9.258-cm diameter horns at FI for a 7.5-deg 
elevation angle. 
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the diffraction from the support tripod, the decision was made to project the tripod blockage onto the 
main reflector similarly to that done with holography. In the area of the projected blockage, the main- 
reector surface currents were set to zero. Using the physical optics scattering program, the fields on the 
subreflector were then computed. Using the reverse scattering program, the focal plane currents in turn 
were computed from the subreflector currents. Using the focal plane analysis, the gain of the antenna with 
support tripod blockage was determined. As a reference, the gain for no blockage was calculated using 
the same technique. For one horn, the loss (or difference in gain) due to blockage was 0.44 dB For an 
array of ™ horns agam the loss due to blockage was 0.44 dB. The array provided an improvement on 
the order of 0.01 dB for both with and without blockage. Therefore, the array was incapable of recovering 
blockage losses, and the calculated losses were due to the support tripod scattering fields outside of the 
region occupied by the array. 


B. Optimization at F3 

The major part of the study was done with the array feed located at F3 since this location would be 
used for a practical implementation. The design approach used followed that used at FI, where first the 
gain as a function of horn diameter was calculated. This was done at one extreme elevation angle, such 
as 7.5 deg, rather than at 45 deg, so as to optimize the array design at the point where the losses would 
e the greatest. This calculation was done at two axial focal positions of 0.0 and —8.89 cm from F3 so 
as to bracket the optimum focal position. Figures 10 and 11 illustrate the results for an array of seven 
horns, where the best horn diameter for a focal position of 0.0 cm was 5.88 cm, and, for a focal position 
of -8.89 cm, the best diameter was 5.08 cm. Next, the calculations were repeated as a function of the 
axial focal position, where the horn diameters were linearly interpolated at each focal position using the 
two best horn diameters previously calculated at the bracketing focal positions. The calculations were 
done at elevation angles of 7.5 and 45 deg. In Figs. 12 and 13, it is seen that the performance for an 
array of seven horns peaks at focal positions of -6.67 cm for an elevation angle of 7.5 deg and -6.1 cm 
for an elevation angle of 45 deg. Since the gain is quite flat in the vicinity of the best focus, a focal 
position between the best focal position for the two elevation angles was selected as the focal position 
for the remainder of the calculations. The selected axial focal position is -6.35 cm for a horn aperture 
diameter of 5.31 cm. 


Having selected the nominal location for a seven-horn array and the best central horn diameter the rest 
of the study consisted of varying the various parameters of the array design. Figure 14 shows the antenna 
performance as a function of elevation angle for the center horn only and for an array of seven horns. At an 
elevation angle of 45 deg, the center horn and the array have the same performance, as expected, since the 
antenna surface was adjusted at this angle. At an elevation angle of 7.5 deg, the single-horn performance 
shows a loss of 1.70 dB relative to the performance at the 45-deg elevation angle. Using the array to 
recover this loss of performance, only a 0.22-dB improvement was obtained. Earlier calculations using the 
orward or transmit mode to compute the possible performance improvement showed improvements of 
about 0.4 dB. The earlier calculations used the array geometry that was used in an experimental program, 
where the array horn diameters were not optimized and the horns were operated 1.67 GHz away from 
their design frequency, allowing more room for improvement. In this study, the horn size was optimized 
as discussed in the previous paragraph. 

The calculation in the first paragraph of this section was for a seven-horn array, where the center 
horn and the six horns in a ring around the center horn had the same diameters. The next calcula- 
tion investigated the effect of using different numbers of horns in the ring around the center horn and 
adjusting the diameter of these horns to completely fill the ring. The larger the number of horns in 
the rings, the smaller the horn diameters. Figure 15 illustrates the performance for this case. In the 
gure, the performance of the center horn is included as a baseline and is a straight line since it is not 
affected by the number of outer-ring horns. Also, the figure shows calculations for one outer ring and two 
outer rings. Where two rings are used, both rings have the same number of horns,’ the second ring nest- 
ing in the first ring. This necessitates the horns in the second ring being larger than those in the first ring. 
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Fi 9* 1 2. Optimum antenna gain versus z (focus) at F3. 
7.5-deg elevation angle (optimum horn diameter). 
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Fig. 13. Optimum antenna gain versus z (focus) at F3. 
45-deg elevation angle (optimum horn diameter). 












For one ring of six horns, the improvement was 0.22 dB, as shown before. However, if the number is 
increased to eight horns, the improvement jumps to 0.63 dB and then drops off for a larger number of 
orns. Using 2 rings, an improvement of 0.72 dB can be gained for 12 horns in each ring. However 
using the second ring, which implies a 25-horn array, only gives an improvement of 0.09 dB over the 
sing e-ring case, which is a 9-horn array. Therefore, the use of a second ring is not practical considering 
the increased complexity. The reason that increasing the number of horns makes an improvement can 
be seen by reviewing the focal plane distribution. Figure 16 shows the focal distribution overlaid with 
the outline of the 5.31-cm diameter horns in the seven-horn case. It can be seen that the center horn 
encompasses the region where the fields are best behaved. The horns in the outer ring, however, are so 
large that they cover an area where the fields slope across the aperture from -5 to -35 dB and therefore 
do a poor job of coupling the fields. The horns like to see a Gaussian distribution for best performance’ 
Figure 17 shows the focal distribution overlaid with the outline of the horns for the eight horns-per-row 
case, where the center horn is 5.31 cm in diameter. Here the smaller horns in the first row do a better job 
of sampling the fields since the fields do not vary more than 10 to 15 dB across their apertures The horns 
in the second row, however, cover regions where the fields are not well behaved and recover very little 
of the field energy. This case shows that increasing the number of horns by two can cause a significant 
improvement. However, it should be noted that the array geometry is driven by the focal distribution, 
which is unique for a given antenna system design and associated aberrations and, therefore the results 
could be significantly different for other antenna designs. This case is interesting because it is ’an example 
of the type of cases that are amenable to the focal plane analysis technique. 


In an earlier study, ray tracing techniques were applied to a seven-horn array, using the same geometry 
used in this article. It was found that with all the array horn axes parallel to each other at F3, at their 
image point at FI, the beams associated with the horns in the ring around the center horn pointed away 
from the antenna axis. This could be likened to an array located at FI having all but its center horn 
rotated outward from the antenna axis. The effect is to improperly illuminate the antenna subreflector. 
At F3, it was found that by rotating all but the center horn inward by 2.62 deg in an aberration-free 
environment, all the beams at FI could be made to be parallel. Another way of viewing this situation 
is that the phase patterns of the focal field distribution at F3 for a BWG antenna are not uniform, but 
tapered. The horns need to be rotated to better match the focal plane phase distribution. Since with 
optimum horn diameters the outer horns are only useful in the presence of antenna aberrations, it was 
of interest to see if rotating the horns could help in recovering more of the losses due to aberrations. At 
an elevation angle of 7.5 deg, a series of calculations was made, rotating all but the center horn inward. 
In Figure 18, it is shown that, for a seven-horn case, rotating the horns 5.2 deg inward improved the 
performance by 0.46 dB. A similar calculation was made for the case that had the optimum eight horns 
m the ring (nine-horn array). In this case, rotating the horns about 4.0 deg had little effect. Evidently, 
significant phase variations in the focal region were beyond the area covered by the smaller outer horns 
for the nine-horn array, but were within the region covered by the larger outer horns of the seven-horn 

array. Figure 19 shows the performance as a function of elevation angle for the seven-horn case with the 
horn rotation angle set at 5.2 deg. 


There are some asymmetries in the BWG geometry that are not significant to horns mounted on the 
antenna optical axis, but could be of concern for off-axis horns such as are used in an array feed system. 
Of concern was what effect these aberrations might have on the array performance as a result of a rotation 
of the antenna in azimuth. Azimuthal rotations cause a change in the relative rotational position of one 
section of BWG mirrors to another. Figure 20 shows the antenna performance at an elevation angle of 
.5 deg as the antenna rotates in azimuth. For the center horn only, the variation is only 0.04 dB, which 
is neg lgible. For the array, the variation is 0.13 dB, still not significant. The other issue is what happens 
when the array is rotated about its axis at F3. Does this affect the way the outer horns sample the 
focal fields? For example, is there an optimum angle? Calculations made on the 10-horn array showed 
variations on the order of 0.1 dB. 
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Fig. 17. The focal distribution overlaid with the outline of the 
horns for the eight horns-per-row case (z (focus) = -6.35 cm, 
elevation angle = 7.5 deg, horn diameter = 5.309 cm). 
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Fig. 20. Antenna gain versus antenna azimuthal angle (z (focus) = -6.35 cm, 
elevation angle = 7.5 deg, horn diameter = 5.309 cm). 

Figures 21 and 22 display the focal plane fields for the antenna and are provided for general interest. 
Figure 21 shows the focal fields for antenna elevation angles of 7.5, 22.5, and 45 deg at the optimum axial 
focal position of -6.35 cm. The primary effect of elevation angle-induced distortions on the focal plane 
distribution is a lowering and spreading of the fields. At 45 deg, the field distribution is approximately 
rotationally symmetric. At FI, the distribution is perfectly circular, as is seen in Fig. 8. The lack of 
perfect circular symmetry shown in Fig. 21(c) is due to the effects of the BWG. Figure 22 shows the effect 
of axial focus changes on the focal fields at an elevation angle of 7.5 deg. There appears to be a rotation 
in the structure of the focal fields in the vicinity of the distribution main lobe as the focus is changed. 
Figure 22(b) shows the best focal position, where the array position is -6.35 cm. 

VII. Analysis of Experimental Configuration 

An experimental program, supported by another task, was performed at DSS 13 using a seven-horn 
array located at F3 that had the same geometry as used in this study, with one exception: The array 
horns were restricted to a nonoptimum diameter of 4.45 cm. Another difference is the dual-mode horns 
used in the experiment were designed for 32.0 GHz but operated at 33.67 GHz. This change in frequency 
degraded the pattern properties of the horns. Because of the degraded horn performance, the results 
presented in the previous section are not typical of what would be expected from the experimental 
program To provide better predictions, different horn-mode models were developed to account for the 
change in performance. While the ratio of the TM„ mode to the TE n mode was 0.405 for the dual-mode 
horns used in this study, to model the horns in the experimental program, a mode ratio of 0.206 with a 
phase of —71.73 deg was required. 

Calculations were made to determine the best axial focal position for the array. Figure 23 shows 
the results for an elevation angle of 7.5 deg. For a single horn, the best position was -7.62 cm, and for 
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Fig. 21. The focal fields for antenna elevation angles of (a) 7.5, (b) 22.5, and (c) 45 deg (z (focus) = -6.35 cm). 

a seven-horn array, the best position was —8.89 cm. Figure 24 shows results for an elevation angle of 
45 deg, and, for both the single horn and for the array, the best focal position is -7.62 cm. Since the 
focus curve is essentially flat in the region of —8.0 cm and the experimental measurements were made at 
-8.89 cm, the predictions were calculated at -8.89 cm. Figure 25 shows the antenna gain as a function 
of elevation angle. At an elevation angle of 45 deg, the gain is 80.66 dB for a single horn and 80.70 dB 
for the array, for an improvement of 0.04 dB. At an elevation angle of 7.5 deg, the gain is 78.89 dB for 
a single horn and 79.39 dB for the array, for an improvement of 0.50 dB as compared with a loss due to 
surface distortions of 1.77 dB. 


To determine if rotations of the antenna in azimuth would affect the performance of the array at F3, a 
series of calculations was made, where the azimuthal position of the antenna was changed in increments 
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Fla 22 The effects of axial focus changes on the focal fields at an elevation angle of 7.5 deg: 
(a) z = 0.0 cm, (b) z = -6.35 cm, and (c) z = -8.89 cm. 


of 45 deg. For an elevation angle of 45 deg, the variation is 0.02 dB for both the center horn and the 
array. The performances for both the center horn and the array are the same since, at an elevation angle 
of 45 deg, the outer horns have very little effect. The variation in performance at an elevation angle of 
7.5 deg is 0.03 dB for the center element and, for the array, 0.13 dB. To determine the best rotational 
alignment of the array at F3, the array was rotated about the optical axis at F3 in fractions of the angle 
subtended by the width of a horn aperture in the ring surrounding the center horn. A series of calculations 
in effect would rotate one horn into the position previously occupied by the next horn in the ring before 
the calculations began. For an elevation angle of 7.5 deg, the gain of the array varied by 0.16 dB. For 
optimal results, the array should be rotationally adjusted. 
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Fig. 25. DSS-13 antenna gain versus elevation angle: focus position, z = -8.890 cm. 


VIII. Conclusions 

The focal plane optimization technique was shown to provide a method of recovering performance 
losses for two antenna problems: (1) the aberration losses due to scanning the beam of an antenna 
and (2) the losses associated with antenna reflector distortions that result from changes in the antenna 
elevation position. The second case really demonstrates the power of the focal plane analysis technique. 
For the antenna configuration associated with the second case, it takes 5.6 h on a Cray Y-MP2 compu er 
to perform one forward scattering calculation to determine the far-field pattern for a single-feed orn. 
The scattering calculation must be repeated seven times to include the effect of each array feed horn for 
a seven-horn array, requiring 39.2 h per antenna configuration. The time to compute the optimum gam 
from the far-field patterns for each feed horn is small compared to the scattering calculation time and will 
be neglected. If the process is repeated for 13 feed sizes, as was done in the example m Section VI. A of this 
article a total time of 509.6 h is needed! The computation time required for the focal plane optimization 
technique can be determined as follows: Assume the time to do the reverse scattering calculation through 
the antenna and beam-waveguide optical system is the same as for the forward scattering case or 5.6 . 

The time required to calculate the currents in the focal plane to the required resolution is 3.3 h. t his 
gives a one-time total time of 8.9 h, which does not have to be repeated unless the antenna geometry 
changes The time required to perform the focal plane optimization that must be repeated for each feed 
size or array geometry ranged from 2 min for the smallest feed size to 33 min for the largest feed size and 

22.5 min for the optimum feed size. The total time to perform the optimization calculation for all 13 feed 
element sizes studied was 2.7 h. The grand total time, including the scattering calculations, comes to 

11.6 h. Reducing the computation time from 509.6 h to 11.6 h illustrates the usefulness of this technique. 
The objection can be raised that you do not need to consider so many element sizes to determine the 
optimum element size. But that is not the issue. There might also be a need to trade off the array 
geometry and/or the horn type, such as single mode, dual mode, hybrid mode, or any other type. 1 is 
could result in the need to analyze more cases than the 13 element sizes considered in the example. 1 he 
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pomt ,s that significant time savings can be obtained by doing the optimization in the focal plane for 
situations where the scattering calculation times are significant, such as is the case for beam- waveguide 

antennas where five or more scattering surfaces must be considered and where their sizes in terms of 
wavelengths are very large. 

The focal plane approach is very accurate. Agreement between the focal plane or reverse scattering 
technique and the classical transmit approach was on the order of 0.02 dB. 
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Significant advances in the development of high electron-mobility field-effect 
transistors (HEMTs) have resulted in cryogenic, low-noise amplifiers (LNAs) 
temperatures are mtUn an order of magnitude of .he no.se 

limit (hf/k). Further advances in HEMT technology at cryogenic temperatures 
may eventually lead to the replacement of maser and superconducting-insulator- 
superconducting front ends in the 1- to 100-GHs frequency band. Key to den to 
fication of the best HEMTs and optimization of cryogenic LNAs are accurate and 
repeatable device measurements at cryogenic temperatures. This article describes 
the design and operation of a cryogenic coplanar waveguide probe system for the 
characterization and modeling of advanced semiconductor transistors at cryogenic 
temperatures. Results on advanced HEMT devices are presented to illustrate the 
utility of the measurement system. 


I. Introduction 

The noise and gain of high electron-mobility field-effect transistor (HEMT) devices have steadily 
improved as the technology is being developed and commercialized worldwide for room temperature 
applications. Although device noise temperatures continue to drop at room temperature, ther^no 
p-narantee that the same improvements will occur at cryogenic temperatures. In order to successfully 
develop ultra-low noise HEMTs for cryogenic applications, one must have quick, accurate, an repea a 
cryogenic data and a good device model. 

The cryogenic on-wafer noise and scattering parameter measurement system ^ rEMT dev^ The 
key to the systematic investigation of the dc and microwave properties of advanced HEMT devices, 
measured parameters are required to evaluate and verify room temperature and cryogenic device models 
for the identification of the optimum cryogenic HEMT structure. The measured parameters along with 
the device models are also needed to design and optimize input, interstage, and output matching circui 
for multiple-stage HEMT low-noise amplifiers (LNAs). 


A cryogenic study and an 
HEMT (PHEMT) device are 


empirical model for the cryogenic operation of an advanced pseudomorphic 
presented. The original apparatus for making these measurements and the 
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II. Cryo-Probe Measurement System 


vanced high-speed transistor technologies, (2) continued advancement of cryogenic LNAs with noTse tern’ 

SsHvbridtnd 6 TtW the qUantUm Hmit {Tn < 5h/A) f ° r grOUnd and ^ace-based applications 
cuits monolithic microwave integrated-circuit (MMIC) semiconductor-superconductor cir- 

“ — * - 

The cryogemc microwave system in this work uses coplanar waveguide probes in a vacuum station 

. g ' & , contains ports for RF cables, thermometers, vacuum pumps, dry nitrogen back-fill lines 
copknar probes with manipulators, and a closed-cycle refrigerator cold beak. The profe bo"y shlwnTn 
g. 1(c), rests on a copper block attached to a fiberglass post. The fiberglass reduces the thermal load 
copper rai ing from the cold head thermally anchors the probe to the 12-K cold station assuring 
sample temperatures of 12 to 20 K. The mechanical and thermal stability of the wafm stl i e^Slshed 

zsr pos,s above the cold “ “ d « — - 

system Tf/7r tant fe f r? ° f thiS design iS the incorporation of a closed-cycle helium refrigeration 

start un Jo * f*?"! ° f s ^ms used open-cycle cooling to reduce 

start-up costs and avoid mechanical vibrations. However, for long term use at the ra fl V , 

down per week a closed-cycle system is significantly less expensive. In fact after the first 6°months 
of operation, the savings in liquid helium is comparable to the price of a tee 

a 3 Im^onarileTots a^ *3°“ "* u° ld “ l ° Probe Sta “°" “» »»mplished with 

two-dimensional bellows and vibration mounts, shown in Fig. 1(b) This svstem «n™„c , lt i ■ , 

microwave measurement from dc 40 GHz over a physical o, S “ 0 K 

he microwave hardware is insulated by vacuum, there am no frosi buildup orTJT .h^tS^T 

M^Wme^um?) aCC ” ate ' “ d “ »— -4. devices 131 

III. Calibration and Measurement Considerations 

(ISS) available from Cascade Microtlh [10 The T RM , USmg an f lm P edance standard substrate 
accuracy than is the SOLT at cryogemc temperatures. In the SOLT method, the short standard 
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introduces uncertainty in the reference plane location because of sensitivity to probe tip placement. The 
LRM method obviates this problem by replacing the short with an open and by having the probe tip held 
approximately 10 mils above the substrate during the calibration sequence. 

The measurement accuracy is also directly related to the calibration conditions. Thermal gradients 
across the gold-plated ceramic probe tips and coax to coplanar transitions [11] alter the electrical charac- 
teristics of the measurement lines. The process of cooling the sample in the laboratory can also produce 
changes in the calibration. The network analyzer typically requires a new calibration if the ambient lab- 
oratory temperature varies by greater than ± 1 deg C. The proximity of the cooled probe system to the 
network analyzer can change the ambient environment (both temperature and humidity). The combined 
results of these effects are appreciable errors in cryogenic temperature measurements. For example, room 
temperature calibrations have been shown to introduce as much as a 20-percent [8] error in the cryogenic 
measurement. In addition, measured results have been reported with room temperature calibrations with 
moding effects [11,12] and deviations from one-pole roll-off [4,11]. 

Since the microwave probe offers a significant thermal load to the device under test, it is evident that 
the chuck and device temperatures are different. For example, with a chuck temperature of 20 K and the 
probes contacting the device under test (DUT), a DUT temperature as high as 50 K has been observed. 
This large temperature differential affects calibration and skews the interpretation of data collected at 
different device temperatures. 


The solution to maintaining calibration integrity and achieving low sample temperatures is to thermally 
anchor the probe body and perform cryogenic calibrations. By thermally anchoring the probe to the cold 
head at 12 K, the thermal load to the DUT is minimized. The remaining microwave hardware [connectors, 
cables, and input to the automatic network analyzer (ANA)] are thermally isolated via vacuum and 
stainless steel hardware. The thermally anchored probe can be calibrated at specific temperatures during 
a measurement cycle. This eliminates the problem of moding (deviation from one-pole roll-off) and allows 
accurate correlation of DUT temperature and measured characteristics. The S-parameter data shown in 
Fig. 2 are an example of data collected in this manner. 

A reliable cryogenic, two-port calibration for S-parameter measurements must satisfy several criteria. 
The return loss of the transmission through standard should be better than -45 dB, and the insertion 
loss should vary within ±0.1 dB of 0 dB. Next, the calibrated open standard, as noted in the previous 
paragraph, should also be within ±0.1 dB of 0 dB, as shown in Fig. 3. A more extensive verification 
can be performed by evaluating additional representative elements not used in the calibration sequence. 
These measurements may include the reflection from a coplanar transmission line (open circuit stub) and 
a 10-dB pad. The input and output scattering parameters of the open circuit stub should be concentric 
spirals on the Smith Chart, with a magnitude of less than 1 and the insertion loss of the pad within ±2 
percent of 10 dB. 

The most accurate and repeatable method of measuring noise parameters at cryogenic temperatures 
is to place the impedance generator within a wavelength of the DUT input. The equivalent noise temper- 
ature of the noise source must also be comparable to the DUT noise temperature. This approach would, 
however, require development of a cryogenic noise generator and noise source. 

For this initial investigation of on-wafer noise parameter measurements at cryogenic temperatures, 
only the probe tips are cooled while the impedance state generator and solid-state noise source (both 
commercially available) are kept at room temperature (several wavelengths away from the DUT). In this 
configuration, the input losses introduce noise comparable to or greater than the noise of the device 
under test (DUT) and reduce the range of available impedance states. For example, in the frequency 
range of 2 to 18 GHz for cryogenic temperatures, the worst-case noise temperature error is ±25 K, while 
device noise temperatures are typically under 10 K. Although this configuration does not provide ac- 
curate single-frequency noise parameter measurements, it does provide for fast and efficient broadband 
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(2- to 18-GHz) on-wafer measurements. The noise calibration must nass all s riQ r. . . 


IV. Small-Signal Modeling 

methods for small signal parameter extraction- m & / pnysics. there are three main 
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Fig. 4. Transmission through a 10-dB attenuator at a physical 
temperature of 20 K using LRM calibration. 



of 3 x 10 12 atoms/cc is grown in a single plane approximately 3 >1 thick. Heavily doped AlGaAs and GaAs 
layers 400 and 350 A thick, respectively, are grown to provide ohmic contacts. A schematic representatmn 
of this heterostructure is shown in Fig. 6. This layered structure produces a conduc ion band discontinuity 
that forms a triangular one-dimensional quantum well at the AlGaAs/InGaAs heterojunction. 
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T-GATE 



Fig. 6. Schematic representation of the PHEMT. 

The noise and scattering parameters were measured for a variety of temperatures and bias conditions 
over a broad frequency span. Scattering parameters were measured from 1 to 25 GHz, while noise 
parameter measurements were attempted from 2 to 18 GHz for the same bias and temperature conditions. 
For all of the experimental conditions, all the S-parameters were measured with an error of less than 2 
percent. A polar plot comparing S 21 at two physical temperatures, 16 and 300 K, is shown in Fig. 7. 
Unfortunately, all four noise parameters could not be reliably measured for the physical temperature and 
frequency spans of interest. At room temperature, all four noise parameters were obtained. A plot from 
2 to 18 GHz showing the scatter in the measured optimum source impedance (real part = RE\Y opt )\ 
imaginary part = IM{T opt \) at room temperature is displayed in Fig. 8. At cryogenic temperatures’ 
however (as discussed in Section III), because of the lossy input and long electrical distance of the 
impedance generator from the device input, only the minimum noise temperature was reliably determined. 
In fact, at the lower frequencies (2 to 8 GHz) and physical temperatures (11 to 39 K) the percentage 
error in the noise parameters was too large (« 80 percent) to be useful. For frequencies above 8 GHz, 
the percentage error in the minimum noise temperature was determined to be less than 25 percent. The 
resulting measurements were then fit to a straight line. The fitting error varied from 6 to 11 percent. 
Figure 9 shows a plot of the measured minimum noise temperature as a function of percentage drain to 
source saturation current at 18 GHz at five physical temperatures, from 16 to 300 K. Over this temperature 
range, the rms error in the minimum noise temperature varied from ±12 K to ± 6 K. 

For this study, the device small-signal circuit elements were determined using two sets of measured 
S-parameters and with the measurement probes calibrated with LRM standards at the physical temper- 
atures of interest. The first set, or cold FET (V d , = 0), measurements determine the bias independent 
extrinsic elements, while the biased device measurements are used to determine the intrinsic elements. 
All of the extracted elements for the different physical temperatures and bias settings are displayed in 
Table 1 along with R t (R s +- R g + R g3 ), estimated T d , and the calculated T min per unit frequency. 

To determine the uniqueness of the extracted intrinsic circuit elements, the extraction algorithm was 
run several times and the results compared. For this test, the extrinsic elements were kept constant and 
initial starting values for a variety of the intrinsic elements were purposely chosen to be up to a factor 
of two larger or smaller than the final value. For all the cases tried, the worst variation of any of the 
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Fig. 9. Measured noise temperature at 18 GHz at several temperatures. 

intrinsic elements was no more than 2 percent. Additionally, S-parameters for a biased device were used 
o simultaneously extract both intrinsic and extrinsic element values. In this case, for the same value of 
fitting error, some of the element values varied as much as 110 percent for R ga or as little as 5 percent 


VI. Analysis and Discussion 

r , P ' 0neering W ° rk ° f Weinreb and Pospieszalski [15,16] has led to the development of HEMT-based 
LWAs that are the lowest noise transistor circuits produced to date. These impressive results have been 
achieved without the benefit of cryogenic, broadband noise measurements. Pospieszalski’s method relies 
on c cryogenic measurements (g m and R da ), manufacturers or measured room-temperature broadband 
S-parameter data (remaining circuit elements), and a single-frequency cryogenic noise-parameter mea- 
surement (i d and T g ) to determine all of the circuit parameters for his empirical noise model. 


For the Pospieszalski noise model [17], the measured noise parameters are required to determine T d , 
e eqtnva ent dram temperature associated with g da , and the equivalent gate temperature T„ associated 
with R gs and as an additional check on the intrinsic element values. Pospieszalski [18] found that for a 
variety of different HEMT and FET devices, T g was equal to the ambient temperature for I ds less than 
20 mA. On the other hand, T d is strongly dependent on I da and, within measurement error, a single 
function may describe its dependence for all devices. (At room temperature, using on-wafer noise and 
scattering parameters from 4 to 18 GHz, Pospieszalski used measured and modeled values for the intrinsic 
ti opt and l min using a least-squares fit to arrive at values of T d and T g .) 


Ideally, for low-noise device studies and circuit modeling, it is preferable to have both broadband 
measured noise and scattering parameters. Since only the measured minimum noise temperatures were 
obtained, only the predicted and measured noise temperatures as a function of drain currents ( I d ,) at 
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Table 1. Extracted elements for the different physical temperatures and bias settings. 


I da a > 

Idas i 

L d , 

Lg, 

L s , 

R d , Rg, Rs , 

Cds > 

Cgd, 

Cg S , 

Gm, 

fid-s > 


fit 

T d , 

Tmin l /, 

% 

mA 

pH 

PH 

pH 

fi n n 

fF 

fF 

fF 

ms 

ft 

ft 

ft 

K 

K/GHz 


T ph y S = 11 K 


15 

2.9 

4 

15 

4 

5.0 

3.8 

1.00 

17.20 

22.90 

76.50 

39.70 

576.70 

1.00 

5.80 

234 

0.123 

25 

4.7 

4 

15 

4 

4.6 

3.8 

1.00 

17.40 

22.20 

84.30 

50.30 

450.00 

1.00 

5.80 

373 

0.153 

50 

9.6 

4 

15 

4 

4.7 

3.5 

1.00 

17.30 

21.60 

91.40 

65.00 

370.00 

1.00 

5.50 

750 

0.196 

75 

14.9 

4 

15 

4 

6.7 

3.0 

1.00 

18.80 

21.20 

92.20 

72.00 

344.20 

0.70 

4.70 

1158 

0.212 

.00 

18.7 

4 

15 

4 

6.0 

3.0 

1.00 

19.00 

19.72 

103.00 

72.00 

331.80 

1.50 

5.50 

1451 

0.292 


Tphys = 16 K 


15 

2.4 

4 

15 

4 

6.0 

3.0 

2.00 

18.80 

3110 

68.60 

36.25 

589.00 

2.00 

7.00 

201 

0.147 

25 

3.9 

4 

15 

4 

6.0 

3.0 

1.00 

17.00 

31.00 

71.00 

45.30 

480.00 

2.00 

6.00 

316 

0.157 

50 

8.4 

4 

15 

4 

6.0 

2.0 

1.00 

17.90 

30.40 

78.70 

60.90 

368.00 

2.00 

5.00 

663 

0.195 

75 

12.7 

4 

15 

4 

6.0 

3.0 

1.00 

15.70 

30.10 

80.80 

68.50 

344.80 

2.00 

6.00 

994 

0.247 

100 

16.4 

4 

15 

4 

6.0 

3.0 

2.30 

17.60 

25.50 

123.60 

67.00 

303.00 

2.00 

7.30 

1279 

0.515 

Tp/iys — 36 K 

15 

2.5 

4 

15 

4 

6.0 

3.0 

2.27 

22.25 

24.77 

72.19 

35.07 

625.75 

11.21 

16.48 

229 

0.381 

25 

4.1 

4 

15 

4 

6.0 

3.0 

2.27 

21.42 

24.41 

80.78 

47.44 

469.76 

7.92 

13.19 

352 

0.403 

50 

8.3 

4 

15 

4 

6.0 

3.0 

2.27 

21.88 

23.39 

92.24 

64.56 

352.77 

6.00 

11.27 

675 

0.500 

75 

12.5 

4 

15 

4 

6.0 

3.0 

2.27 

15.47 

29.67 

90.21 

74.60 

298.41 

5.50 

10.77 

999 

0.547 

100 

16.4 

4 

15 

4 

6.0 

3.0 

2.27 

16.91 

24.56 

133.61 

66.91 

303.99 

5.50 

10.77 

1299 

1.021 


Tphys = 49 K 


15 

2.2 

4 

15 

4 

6.0 

3.0 

2.27 

21.83 

24.93 

72.19 

35.42 

526.80 

10.88 

16.15 

218 

0.464 

25 

3.7 

4 

15 

4 

6.0 

3.0 

2.27 

21.51 

24.39 

71.92 

47.40 

468.31 

7.20 

12.47 

334 

0.398 

50 

7.5 

4 

15 

4 

6.0 

3.0 

2.27 

21.88 

23.39 

92.24 

64.56 

352.77 

6.00 

11.27 

627 

0.562 

75 

11.1 

4 

15 

4 

6.0 

3.0 

2.27 

20.65 

23.06 

96.51 

72.77 

324.14 

6.00 

11.27 

904 

0.654 

100 

15.8 

4 

15 

4 

6.0 

3.0 

2.27 

21.69 

17.45 

137.63 

66.24 

315.41 

6.00 

11.27 

1266 

1.229 


Tphys = 300 K 


15 

2.7 

4 

15 

4 

6.0 

3.1 

3.40 

17.10 

19.50 

69.60 

28.13 

798.00 

2.00 

8.50 

508 

1.253 

25 

4.6 

4 

15 

4 

6.0 

3.0 

4.90 

19.30 

18.33 

83.60 

40.36 

560.00 

2.00 

9.90 

654 

1.533 

50 

9.1 

4 

15 

4 

6.0 

1.5 

3.40 

18.70 

18.02 

90.00 

51.00 

464.90 

2.00 

6.90 

1001 

1.480 

75 

13.8 

4 

15 

4 

6.0 

3.0 

3.50 

18.40 

17.41 

96.60 

57.90 

417.70 

2.00 

8.50 

1363 

1.912 

100 

18.4 

4 

15 

4 

6.0 

1.96 

3.70 

18.24 

17.05 

101.60 

60.40 

398.30 

2.00 

7.66 

1717 

2.104 
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18 GHz for a variety of different temperatures are compared, 
is calculated using Pospieszalski’s expression, 


The predicted minimum noise temperature 


-i _ ^7T fCg S yj ^JdsR-tTgTd 
min — 

9m 


for the minimum noise temperature. The small-signal circuit elements are those extracted from the 

measured S-parameters, while T g is the device physical temperature and T d is estimated from the results 
in [loj. 


Figure 10 shows a comparison of the room temperature results as a function of percentage of drain 
saturation current. The measured minimum noise temperature is more than a factor of two larger than the 
pre icte va ues. This discrepancy is probably associated with an unrealistic extracted value for R t . The 
extracted value for R t at room temperature is approximately the same as the value at 11 and 16 K and 
50 percent smaller than the 36- and 49-K extract value (see Table 1 and Fig. 11). The room temperature 
^’.should have been 2 to 4 times the cryogenic value. However, at cryogenic temperatures for low 
ram-current values, the measured and modeled results agree within the measurement uncertainties (see 


For single-frequency measurements at 18 GHz, the worst-case errors for the measured minimum noise 
temperature are ± 25 K and ± 15 K for the calculated model values. That the cryogenic results agree 
is possibly fortuitous. However, that there is agreement indicates that a closer examination with a 
more accurate and reliable method of cryogenic noise parameter measurement and small-signal element 
extraction technique is absolutely necessary. 


Three of the four intrinsic circuit elements, g m , g ds . and C g ,„ that determine T mtn behaved as expected, 
e transconductance, g m , and the output conductance, g ds , were observed to show significant increases 
upon cooling and increasing drain current (see Figs. 13 and 14). The C gs shown in Fig. 15 is relatively 
constant with temperature but increases linearly with increasing bias while, on the other hand, the gate 
resistance, R gs , flip-flops, increasing in resistance as the temperature is lowered toward 36 K but'dropping 
o its room temperature value at 16 K and half that value at 11 K (see Fig. 11). In addition, for high 
rain-current values, the resistance appears to be independent of drain current value Clearly more 
accurate noise-parameter measurements at these cryogenic temperatures will significantly increase our 
HEMTs 56 ° f the temperatUre behavlor of the parameters that determine cryogenic noise performance of 


VII. Conclusion 

The feasibility of cryogenic, broadband on- wafer scattering parameter measurements has been demon- 
strated. It was also shown that a commercial noise- parameter measurement system is inadequate for 
cryogemc applications. Although some useful cryogenic noise temperature data were obtained, it was not 
o sufficient density or accuracy for a systematic device study. Systematic studies are necessary to develop 
device technologies for future applications and to properly select, develop, and design high-performance 
cryogenic LNAs. 


The development of a cryogenic probe tip integrated with a cryogenic noise and impedance generator 
capable of performing both scattering and noise parameter measurements will circumvent the limitations 
posed by the current characterization techniques. This type of probe will allow nondestructive char- 
acterization of devices and circuits while still in wafer form. This measurement technique will provide 
signfficant cost savings by improving accuracy and design methodology and by reducing test verification 
c sts. Once operational, such a system will be beneficial in studying HEMT LNAs as well as hybrid 
semiconductor/superconductor circuits. 
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l P UI?K , T . MT d and LNAs with maser-like noise temperatures. The advent of suDer 

mullT , HEM y t 7“ S Wi " “ t0 g ' eater f "*W coverage and economical cryogej^ stng^Zd 
multiple-element feed systems for routine telemetry and radio navigation. ~ 
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Radio-Science Performance Analysis Software 

D. D. Morabito and S. W. Asmar 
Telecommunications Systems Section 


The Radio Science Systems Group (RSSG) provides various support functions 
for several flight project radio-science teams. Among these support functions are 
uplink and sequence planning , real-time operations monitoring and support, data 
validation, archiving and distribution functions , and data processing and analy- 
sis. This article describes the support functions that encompass radio-science data 
performance analysis . The primary tool used by the RSSG to fulfill this support 
function is the STBLTY program set. STBLTY is used to reconstruct observable 
frequencies and calculate model frequencies , frequency residuais, frequency stabil- 
ity in terms of Allan deviation, reconstructed phase, frequency and phase power 
spectral density, and frequency drift rates. In the case of one-way data, using an 
ultrastable oscillator (USO) as a frequency reference, the program set computes the 
spacecraft transmitted frequency and maintains a database containing the in-flight 
history of the USO measurements. The program set also produces graphical dis- 
plays. Some examples and discussion on operating the program set on Galileo and 
Ulysses data will be presented. 


I. Introduction 

The radio science investigations of the different flight projects require a diverse set of tools for ensuring 
that acquired radio science data satisfy experimenter requirements. For example, the gravitational wave 
experiments on Galileo, Ulysses, Mars Observer, and Cassini each has specific stability requirements on 
the received signal frequency and amplitude. The Radio Science Systems Group (RSSG) utilizes various 
software tools to verify that the received data satisfy the requirements and that the required configurations 
are in place. STBLTY is the primary analysis tool used by the RSSG to characterize the performance 
of the radio science instrument. The STBLTY program set has been utilized by the RSSG to perform 
measurements of ultrastable oscillator (USO) frequency for the Galileo Solar Redshift Experiment [1,2], to 
characterize the performance of the Galileo Orbiter USO [3], to verify performance of the Mars Observer 
USO, 1 and to verify performance for the gravitational wave experiments on Galileo [2] and Ulysses [4]. 
The characterization of the Galileo USO is important for occultation radio science experiments that will 
be performed during the Jupiter tour [5]. This article will focus on describing the algorithms that are 
performed by the individual component programs of the STBLTY program set and will present examples 
using STBLTY on Galileo and Ulysses radio science data. 

The version of STBLTY used for Voyager radio science was inherited by the Galileo and Ulysses 
radio science projects. This version was used primarily for one-way data occultation analysis, where 

! D. D. Morabito, “Mars Observer 93-049 USO Test Processed Through STBLTY;’ JPL Interoffice Memorandum 
3394-93-077 (internal document), Jet Propulsion Laboratory, Pasadena, California, June 17, 1993. 
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a baseline frequency signature was used to remove all effects other than those due to the planetary 
atmosphere. Hence, the inherited version did not have the stringent accuracy requirements needed for 
Galileo and Ulysses radio science experiments. In order to make STBLTY suitable for Galileo and 
Ulysses cruise data, several upgrades were made to the program set after the Voyager era (Neptune 
1989). These upgrades included adding the capability of computing coherent two-way Doppler residuals, 
modeling the effects of a spinning spacecraft, upgrading the troposphere calibration model, introducing 
an ionosphere calibration algorithm using tracking system analytic calibration (TSAC) data as input, and 
estimating a more accurate trajectory model. The new trajectory model includes effects due to nutation, 
irregularities in Earth rotation (UT1-UTC), polar motion, equation of the equinoxes, zonal tides, and 
more comprehensive formulation of state vectors. 

Figure 1 is a block diagram displaying the interconnections between the component programs of the 
STBLTY program set. Input radio science data consist of closed- loop data from the DSN tracking system 
in the form of archival tracking data files (ATDFs) or open-loop data from the radio science system in 
the form of original data record (ODR) tapes. The input celestial reference set (CRS) file from the 
Navigation Team of a supported project provides the spacecraft-state vectors necessary to model the 
spacecraft trajectory from the observable data. 


INPUTS 


OUTPUTS 



Fig. 1. The STBLTY program set block diagram displaying program names, Input, output, 
and intermediate tapeffile products. 


The program TRAJVECT reads in body-centered spacecraft state vectors from the CRS file, translates 
the Earth-centered vectors to the location of the observing station, and performs light-time solutions. 
TRAJVECT produces a set of downlink state vectors at the time of signal reception (in file F45), and if 
data are two-way, produces a set of uplink vectors at the time of signal transmission (in file F49). These 
data sets are passed to the program RESIDUAL. 

The closed-loop tracking system utilizes a phase-locked loop (PLL) that performs signal acquisition, 
lock-up, and detection in real time. The closed-loop data are conditioned by the Multi-Mission Navigation 
Team, which produces and delivers an ATDF to the RSSG. Among the quantities on an ATDF are Doppler 
counts, Doppler “pseudoresiduals” (residuals based on predicted frequencies used to tune the receivers), 
signal strengths (AGCs), and Doppler reference frequencies either in the form of a constant frequency 
or uplink ramps. OBSERVE reads the Doppler extractor cycle counts and Doppler extractor reference 
information from the ATDF and reconstructs these into received sky frequencies (in file F50). If data are 
two-way, OBSERVE writes the uplink ramps into the programmable ramp file (PRF). 
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The open-loop system mixes an incoming intermediate frequency (IF) signal with a signal whose 
frequency is a linear approximation of the predicted frequency, also known as the programmable oscillator 
control assembly (POCA) ramps. The baseband signal is passed through a filter whose bandpass is 
centered at the expected frequency and has a sufficiently wide bandwidth to allow for any unexpected 
signal frequency excursions. A set of analog-to-digital converters (ADCs) digitizes the received bandwidth 
and then writes the samples onto the ODR tapes. The program PLLDEC reads the ODR tapes and 
performs signal detection on the recorded samples, outputs an F36 file with the detected frequencies and 
signal strengths, and extracts the open-loop receiver tuning information (POCAs) from the ODR headers 
and writes these to an F33 file. The program OBSERVE reads the PLLDEC F36 and F33 output files 
and reconstructs the information into received sky frequencies and writes these in an F50 file that is 
passed to RESIDUAL. 


The program RESIDUAL reads the observed sky frequency F50 file and, if the data are two-way, the 
PRF file output from OBSERVE and the spacecraft trajectory vectors from the TRAJVECT F45 (one- 
way or two-way) and F49 (two-way only) output files. RESIDUAL computes model frequencies from the 
F45 and F49 trajectory vectors and a spacecraft spin model. These model sky frequencies are subtracted 
from the observed sky frequencies to produce residual frequencies that are written onto the RESIDUAL 
output F52 file. 

The F52 residual frequencies are input to the STBLTY program, which performs further corrections 
and produces Allan deviations, frequency and phase power spectral densities, and other quantities. In 
the case of one-way USO data, STBLTY computes the spacecraft transmitted frequency and writes out 
records for each pass onto a database summary file, SUMFIL, which is accessed by other programs for 
USO measurement analysis. The SUMFIL is also delivered to the radio science investigator, who performs 
an analysis for the purpose of measuring the solar gravitational redshift [1], 


II. Program Descriptions 

A more detailed description of each of the component programs of STBLTY follows. 

A. PLLDEC — Open-Loop Data Signal Detection 

The scheme employed for detecting signal frequency and power from the digitized samples of the 
open-loop system is a digital second-order phase-locked loop program called PLLDEC. The digitized 
samples output from the ADCs in the Deep Space Communications Complex (DSCC) spectrum processing 
(DSP) subsystem are written onto the ODR tapes (see Fig. 1) that are input to PLLDEC. A theoretical 
background that supplements the following discussion is given in [6], and a detailed description of a 
similar version of PLLDEC (using noncoherent AGC requiring high signal-to-noise ratio (SNR) signals) 
is contained in an internal memorandum. 2 

PLLDEC reads the set of digitized samples, x(i), i = !,■■-, N, from the ODR tape. An initial fre- 
quency, /o, for the PLL is obtained from a fast Fourier transform (FFT) performed on an initial subset 
of samples. The output signal frequency, y(nT), is expressed as a linear combination (recursive function) 
of the recorded samples, x(i), and the previous output signal-frequency estimates, y(iT). 

The implemented phase-locked loop consists of a mixer (phase detector) and loop filter providing the 
forward gain path, F(s ) = (1 + as)/bs; a low-pass filter, K(s), for suppressing the mixer sum product 
(with negligible baseband response); and an integrator, G(s) = 1/s. The closed-loop transfer function 
(the transform of the PLL phase over the transform of the input phase) is given by 


2 A. Densmore, “A Digitally Implemented Phase-Locked Loop Detection Scheme for Analysis of Phase and Power Stability of 
a Calibration Tone,” JPL Interoffice Memorandum Voyager RSST-88-016 (internal document), Jet Propulsion Laboratory, 
Pasadena, California, February 4, 1988. 
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rM F(s)G(s) 1+as 

[S) 1 4- F(s)G(s) 1 + as 4* 6s 2 


The desired s-domain transfer function of the second-order PLL of the output detected frequency to the 
input phase is given by 


W(s) 


s(l 4- as) 

1 -h as 4- bs 2 


The continuous time analysis of the discrete time PLL is justified by assuming that the product of the 
PLL noise bandwidth and sample period is much less than unity, (B PLL T « 1). In order to realize the 
digital PLL implementation, the bilinear transformation from s-space to z-transform space is employed. 


It has been shown that this transformation outperforms others for the case of low sample rates [7]. 

By substituting Eq. (1) for s into G{s), the following z-domain expression for the integrator is obtained: 


G(z) = 


T 

2 



( 2 ) 


By substituting Eq. (1) for s into F{s), the corresponding z-domain function for the loop filter is obtained 

T + 2a [1 + - 2a)/(T + 2a))] ^ 

F ^ Z > " 2 b 1 -z- 1 


The relations for a and b are given by 


R+ 1 

4 BpLL 



where R is the loop damping parameter, which is set equal to 2, the optimum value of loop performance 
specified by Jaffee and Rechtin [8] when the initial phase is unknown but uniformly distributed. A 
typical value for the one-sided loop bandwidth used for Galileo and Ulysses radio-science data processing 

is Bpll — 1 Hz. 


A single-pole low-pass filter is also employed to suppress the mixer sum product and has the following 
s-domain representation: 


K(s) = 


1 

1 4- ( s / a) 
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Using the bilinear transformation, the implemented z-domain transfer function of the low-pass filter is 
given by 


K(z) 


aT 1 + z -1 

2 + aT [1-z-i ((2 - aT )/( 2 + aT))} 


( 4 ) 


where a is set equal to 10 B PLL . The bilinear form of Eqs. (2), (3), and (4) are used to determine the 
appropriate coefficients that multiply the previous input samples and previous output values to yield each 
output estimate of signal phase or frequency. 

The contribution to the Allan deviation due to additive white Gaussian noise (AWGN) over the PLL 
bandwidth B PLL can be estimated from the following equation: 




v/3flp LL CiVfl-i 

27t/ 0 t 


where -Bpll is the PLL bandwidth, Hz; CNR is the carrier signal-to-noise ratio of the signal; /o is 
the observing frequency, Hz; and r is the time interval, s. By knowing the expected signal strengths 
and system noise temperatures, one can use this equation to estimate the B PPP required to obtain an 
estimate of the AWGN contribution to the Allan deviation. By reducing B PLL appropriately, one can 
reduce this contribution below other effects, thereby allowing visibility of those effects. Typical CNRs for 
Galileo and Ulysses result in AWGN contributions to the Allan deviations at 1000 s, which lie well below 
the measured values that are limited by media (for two-way Galileo and Ulysses) and USO (for Galileo 
one-way) when using 1 Hz for B P n. 

B. TRAJVECT— Trajectory Vector Conditioning 

The program TRAJVECT reads spacecraft state vectors from the CRS file, corrects Earth-center 
vectors to the observing-station location, and performs light-time solutions to yield a set of downlink 
state vectors at the time of signal reception, and (if two-way) performs a light-time solution to yield a 
set of uplink state vectors at the time of signal transmission. The conditioned set of state vectors is then 
passed to RESIDUAL. For discussions on how the NAV orbit determination solutions are performed, 
from which the CRS files are derived, the reader is referred to [9] for Galileo and [10] for Ulysses. 

1. Input CRS File State Vectors. For cruise processing, where the gravitational effects of the 
planets are assumed negligible, it suffices that the CRS file contains EME 1950 (Earth mean equator and 
equinox of B1950.0) state vectors of the spacecraft relative to the Sun and the Earth, as follows: 

x ®- 8 /c(<fc) = Earth center-to-spacecraft position vector 
v ® s/c {tk) _ center-to-spacecraft velocity vector 

x O-s/c ^ k ) _ Sun-to-spacecraft position vector 
v © s/c (tk) _ Sun-to-spacecraft velocity vector 

where k is the index of the vectors which are evenly spaced in time, T s apart. The time tag, t k , is 
referenced at the spacecraft in ephemeris time in seconds past 1950.0. 

2. Conversion to Sun Center. The state vectors from the CRS file are referred to Sun center using 
the following translations: 
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x®-®(tfc) = x®~ s/c (t k ) - x®-* /c {t k ) (5) 

x s/c - ®(t fc )= - x°" s/c (tfc) (6) 

v®-®(t fc ) = V®~ a/C (t k ) (j-^j - v©- s/c (t fc ) (7) 

v s/c_0 (<fc) = -v 0 - 9 / c (t fc ) (8) 

where the (+/rel) notation means that the velocities add relativistically. 


3. Calculation of Observing Station Vectors. TRAJVECT translates the Earth-Sun vectors 
[Eqs. (5) and (7)] to the location of the observing station. The station-Earth vectors are computed in 
the EME 1950 frame and are added to the Earth-Sun vectors [Eqs. (5) and (7)] to produce the required 
station-to-Sun state vectors. Among the corrections applied to the station vectors are precession, nutation, 
UT1-UTC, polar motion, equation of the equinoxes, and solid Earth tides. 

A table of station locations in cylindrical coordinates (referenced to the mean pole, equator, and prime 
meridian of 1900.05) are read to obtain the observing station coordinates: radius off the spin axis, r s 
(km); height above the equatorial plane, z (km); and longitude in degrees east of the Greenwich meridian, 
A. Thus, the body-fixed station vector in rectangular coordinates is 

r s cos A 
x dss = r 3 sin A 

z 

Several notations are used to denote time: the corresponding Julian date of t k , 

* 86T400 + 2,433,282,5 

the number of Julian centuries of 36,525 mean solar days that have elapsed since 1950.0, T k , 

k 3,155,760,000 

and the difference of ephemeris time (ET) and universal time (UT1), A t k . 

To translate the station vectors to the frame of the current epoch, several translations and rotations 
are applied. Polar motion matrix rotations are applied, and the vectors are then rotated to the correct 
position in the frame of the current epoch. The rotation matrices that perform these corrections use values 
of universal time and polar motion obtained from the monthly publications of the International Earth 
Rotation Service (IERS). These values are interpolated to the appropriate time tag, t k , of the vectors, 
yielding the Earth orientation values of AturiR (UT1R-UTC) and polar motion X p and Y p . The AturiR 
is further corrected by adding the effects of short-period zonal tides (monthly and fortnightly), to yield 
Aturi- The zonal tides are computed using the formulation of Yoder et al. [11]. 

The angle of the Greenwich meridian at the current epoch, corrected for variation in Earth rotation, 
<*g( 4)> is known as Greenwich mean sidereal time, or the Greenwich hour angle of the mean equinox of 
date, and is given by 
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(*G(tk) — ao(tk) + cti(tk) + a 2 (tk) + a3(tk) 


(9) 


where a 0 (t k ) = -(1.002737909294 + 7*0.589 x 10~ 10 )((Af*)/240) + 100.0755426042 + T fc 2 0.0003870833 is 
the fraction of a turn (in degrees) of the Earth over 7* Julian centuries, c*i(f*) = MOD[(tk/ 240), 360) is 
the fraction of a turn (in degrees) due to daily rotation, a 2 (t k ) = MO£>[36, 0007*, 360] is the fraction of 
a turn (in degrees) due to the revolution of the Earth in its orbit, and q 3 (i*) = 0.76931208337* is the 
fraction of a turn (in degrees) over T k Julian centuries, which accounts for the slowing down of the Earth 
due to the exertion of tidal forces. 

The nutation model used by the program set is the International Astronomical Union (IAU) 1980, 
or Wahr nutation model, which was adopted by the IAU in 1984 and uses a series developed by Wahr 
based on the nutation model of P. K. Seidelmann [12]. The expected accuracy of this model (1 milliarc) 
[13], translates to about a 1-mHz error at 8.4 GHz, but is expected to be constant over a tracking pass. 
The angles of nutation estimated from this model are Ae, the nutation in obliquity; A’i', the nutation in 
longitude, and e m , the mean obliquity of the ecliptic. The derivatives are computed from the differences of 
the nutation angles evaluated at two slightly different times divided by the time difference. The obliquity 
is added to the mean obliquity to yield the true obliquity of the ecliptic, e. 

The Greenwich mean sidereal time in Eq. (9), a G (t), is converted to Greenwich apparent sidereal time, 
a GA (t) (the Greenwich hour angle of the true equinox of date), by applying the correction of the equation 
of the equinoxes: 


&GA(tk) = OLc(tk) + A# cos e 


The space-fixed station vector projected in the celestial frame at t k (the position of date) is computed 
from the body-fixed station vector by applying a series of rotations that corrects for polar motion (Rj and 
R 2 ) and Earth rotation (R 3 ). This ' true of-date” space-fixed station vector is corrected for nutation (N), 
which performs the rotation from the true of-date equator and equinox to the “mean of-date” equator and 
equinox. Finally, the station vector is rotated from the mean of-date equator to the mean equator of the 
reference epoch (EME 1950) using the precession rotation matrix, P. See the Appendix for definitions of 
P and N. The rotation matrices Ri,R 2 , and R 3 are given by 



'cos a G A(tk) 

— sin < 

R3(«G/»(^)) = 

sin a CA (t k ) 
0 

COS Ol 


f cos X„ 0 

— sin X. 

R 2 (A p ) = 

0 

1 

0 


L sin X p 0 

cos Xp 


'l 

0 

0 

Ri (Y P ) = 

0 cos Y p 

sin Y p 


0 - sin Y p 

cos Y p 


0 

0 

1 


The station position vector in EME 1950 coordinates is computed as follows: 


x ds9 “ e (t fc ) = PN[R 3 R 2 R 1 x d8S + Ax, dcs ] 


( 10 ) 
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where R x , R 2 , and R 3 are the “Earth platform” rotation matrices given above, N is the nutation rotation 
matrix, P is the precession rotation matrix, and A x. tides (t k ) designates the perturbations due to solid 
Earth tides which, along with velocity perturbations, Av ti d es , are induced by the Sun and Moon. 

The station velocity vector is given by 4,5,6 

v dss "®(f fc ) = [PN + PN][R 3 R 2 Rix dss + A x t ide»(tk)] + PN[R 3 R 2 Rix dss + A v Hdes (t k )] (11) 

where the N and P matrices (see the Appendix) denote the derivatives of the nutation and precession 
matrices, respectively, A v tides (t k ) is the station velocity perturbations due to tides, and R 3 is given by 


R3(<*GA(ffc)) 


— sin ctGA(tk) —cos atGA(Tk) 0 
cos acA^k) -sin a G >i(f*:) 0 

0 0 0 


4. Light-Time Solutions. The light-time solution is performed as follows. First, an initial estimate 
of the one-way light time is made using Eqs. (5) and (6): 

, |x®-°(f*)-x 9 /<=-®(t fc )| 

r 0 (tjfc) = — 


where the “0” subscript denotes the initial estimate. The ith estimate of the time tag of downlink 
reception a one-way light time later is made as follows. 

tk,i = t It + Ti(tk,i- 1) {tk,0 = ffe) (12) 


The three sets of vectors with time tags, and £/+i, lying closest to the estimated received time, 

t ki , are quadratically interpolated to yield estimates of the state vectors at t kji . The station-Earth 
vectors [Eqs. (10) and (11)] are evaluated at t k ,i and added to the Earth-Sun vectors (Eqs. (5) and (7) 
interpolated at tk,i) as follows: 


x dss-0 


v dss-0 


(fjfc,i) = X ds * - ®(ffc,i) +X® - ®(tfc,,) 

(13) 

(t kti ) = v dss- ®(ffc,i) (dj) v®-®(t fcli ) 

(14) 


The correction to the light-time of the next iteration is 


3 T D Moyer “Correction to Earth Fixed Station Coordinates Due to Solid Earth Tides, Ocean Loading, and Pole Tides 
and Calculation of Periodic Terms of UT1,” JPL Engineering Memorandum 314-505 (internal document), Jet Propulsion 
Laboratory, Pasadena, California, July 4, 1991. 

4 N. A. Mottinger and T. D Moyer, “A Close Encounter With a Doppler Observable,” JPL Interoffice Memorandum 
314.5-1025 (internal document), Jet Propulsion Laboratory, Pasadena, California, June 18, 1984. 

5 t D Moyer “Proposed Changes to ODP Transformation Between Body-Fixed and Space-Fixed Coordinates for the 
Planets and the Sun,” JPL Engineering Memorandum 314-271 (internal document). Jet Propulsion Laboratory, Pasadena, 
California, June 16, 1982. 

6 T. D. Moyer, “Time Derivative and Partial Derivatives of the Body-Fixed to Space-Fixed Coordinate Transformation for 
the Sun, Planets and Planetary Satellites,” JPL Engineering Memorandum 314-379 (internal document), Jet Propulsion 
Laboratory, Pasadena, California, August 22, 1985. 
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r i + l(£fc,i) — r i(tk t i-l) + Ar l+ i(^,i) 


The time tag for the next iteration is refined as [Eq. (12), substituting i + 1 for i] 


— tk + 'Ti+l(^ t i) 


The correction [Eq. (15)] is compared against a preset tolerance that should be tight enough for 
accurate vector determination given the stability requirements of the experiment being processed and the 
time spacing of the input CRS vectors. The current value used for Galileo and Ulysses cruise data is 
10 s, assuming vectors separated by 30 s. Once the criterion is met, the iterative solution is completed 
and returns the estimated vectors at If the criterion is not met, then additional iterations are 

performed until it is satisfied. In the case of the uplink vectors, the procedure is the same except that a 
minus sign replaces the plus sign in the appropriate equations. 

The time tag of the final iteration of the light-time solution for the downlink output vectors is redefined 
for notational convenience as 


*k ~ = £fc,i + l 


The M state vectors from the CRS file are processed this way such that a set of M heliocentric station 
and spacecraft state vectors referenced at times k = 1, ■ ■ • , M is produced in the EME 1950 reference 
frame. These vectors are passed to RESIDUAL via the downlink vector file F45 and uplink vector file 
F49. For downlink, 

*dn~ G W) 

v£T“ g WE) 

*Xr 0 «> 

v£-®«) 

For uplink, 

x up“°(£f) = station position 
v up S-G (*f) = station velocity 
x up C G (tj*) = spacecraft position 
Vup C °(i*) = spacecraft velocity 

where the subscript k is the index of the kth vector in the F45 downlink vector file, and the subscript l 
is the index of the Zth vector in the F49 uplink vector file. 


= station position 
= station velocity 
= spacecraft position 
= spacecraft velocity 


129 



C. OBSERVE— Sky Frequency Reconstruction 

1. Closed- Loop System. The program OBSERVE reads in “Doppler” cycle counts from the closed- 
loop system, <p(tj), along with information needed to compute the Doppler extractor reference frequency, 
fref{tj)- The cycle counts consist of an integer number of counts output from the Doppler counter, plus 
a fractional term output from the Doppler resolver. The Doppler reference frequency may either be a 
constant synthesizer frequency (SIM), normally used as the reference during one-way data acquisition, 
or be computed from uplink ramps used as the reference during two-way data acquisition. The ramps 
used to tune the uplink frequency consist of a series of programmable frequencies and ramps that drive 
a digitally controlled oscillator ( D C O) , which is multiplied up and input to the exciter that drives the 
transmitter. The ramps are written to the PRF file, which is input to RESIDUAL, where they are used 
to reconstruct the uplink frequencies for two-way data. 

For closed-loop data read off an ATDF for a standard Deep Space Station (DSS), OBSERVE computes 
the received biased Doppler frequencies from the accumulated cycle counts as follows: 


fdop(t. j) 


+ 1) 0fe-l) 

tj + 1 tj— i 


In the case of a modulo reset, where <f>{tj) is less than 1 ), 


. , <t>(t i+ i) + 2 32 

JdovV'j) — . , 

Cj + i tj - 1 

The Doppler frequency is converted to a 2.3-GHz sky frequency by 

= 96 22 ^fref(tj) —B ( fdop(tj ) — 10 ) 

where 10 6 is the 1-MHz bias, B is the sign of the bias, f ref is the Doppler extractor reference frequency 
interpolated to the time tag of the measured Doppler, and 240/221 is the spacecraft transponder ratio. 
For 8.4 GHz, 


riW = 96 —B (/**&) - 10 6 ) 

In the case of the 34-m high-efficiency (HEF) station, the 8.4-GHz frequency is given by 7 

fsky{tj ) = [32 (4.68125/re /(f; ) — 81.4125 10 6 ) + 6500 10 6 ] —B ( fdop(tj ) — 10 6 ) 

2. Open-Loop System. For the case of open-loop data, the frequencies detected from the recorded 
samples by PLLDEC, f rec , in the F36 file are converted to sky frequencies, f 3 or / x , using the tuning 
frequencies, f poca , interpolated from the POCA offsets and ramps in the F33 file at the f rec time tag. 
For the radio science IF-VF converter assembly (RIV) receivers, 


7 T. D. Moyer, “Change to the ODP and ODE for Processing X-Band Uplink Data,” JPL Engineering Memorandum 314-430 
(internal document), Jet Propulsion Laboratory, Pasadena, California, October 15, 1987. 
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f s (tj) = 3 (fpocaitj) + (™ X 10 6 ^ + (1950 x 10 6 ) - f 3amp + f rec ( tj ) 

f (tj) — 11 ifpoca(tj) — (10 X 10 6 ) 4- (8050 X 10 6 ) — 3 f 9amp + frec(tj) 
where f SQmp is the sampling frequency of the recorded samples. For the multimission receivers (MMRs), 


/ i^j) — 48 fpoca(tj) + (300 X 10 6 ) — -fsamp + frec(tj) 



48 fpocaitj) + (300 x 10 6 ) 


3 

4 


4 " frec{tj) 


The observed sky frequencies from OBSERVE are input to RESIDUAL, where model and residual fre- 
quencies are calculated. 

D. RESIDUAL — Residual Frequency Calculation 

The program RESIDUAL reads the sky frequencies from the OBSERVE output F50 file and the 
state vectors from the TRAJVECT output F45 and F49 files to compute model frequencies based on 
the spacecraft trajectory, and produces frequency residuals for either 2.3 GHz, 8.4 GHz, or differenced 
S-3/11X data. RESIDUAL applies a polarization-dependent correction for a spinning spacecraft 8 and 
can estimate and remove a sinusoidal signature induced by an off-axis antenna on a spinning spacecraft 
(e.g., low-gain antenna (LGA)-2 on Galileo). 

1. Spin Correction. For the constant spin correction, specific default values are built into the code 
to account for different spacecraft configurations. This correction assumes the 2.3-GHz one-way default 
values for f spin shown in Table 1 (note that Galileo has a 3-rpm spin rate and Ulysses has a 5-rpm 
spin rate). This correction could also be obtained from attitude-control telemetry data and the defaults 
therefore overridden. For Galileo passes involving LGA-2, the spin rate could also be solved for from a 
three-parameter sinusoid fit of the trajectory corrected residuals. 

Table 1. One-way default values for f spin . 


Spacecraft 

Antenna 

Spin mode 

ff P ln WaV < Hz 

Galileo 

LGA1, HGA 

Dual 

0.052 


LGA2 

Dual 

-0.052 


LGA1, HGA 

All 

0.048 


LGA2 

All 

-0.048 

Ulysses 



-0.0833 

Any other 



0.00 


RESIDUAL applies the appropriate spin correction to the sky frequency depending upon the spacecraft 
spin configuration and frequency band. If the data are 2.3-GHz one-way, f spin is applied as given in 
Table 1. For the case of two-way data with a 2.3-GHz downlink and a 2.3-GHz uplink, 


8 P. Priest and J. Breidenthal, “Galileo Spin Effects,” JPL Interoffice Memorandum 3393-90-186 (internal document), Jet 
Propulsion Laboratory, Pasadena, California, November 30, 1990. 
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rS/S _ f 240 i_ lA fSA-way 
J spin ~ l 221 ^ 1 / spin 

For the case of two-way data with an 8.4-GHz downlink and a 2.3-GHz uplink, 


fS/X 

J spin 


240 11 „ 

— 1 

221 3 


fS,l — way 
f spin 


For the case of dual-frequency S-X residuals, the spin correction applied is given by 

rS-X _ rS, 1-way 
J spin 11^ s P* n 

2. One-Way Data Processing. The state vectors from the TRAJVECT F45 file are interpolated 
to the time tag of the observable frequency (tf ss at the observing station and tj at the spacecraft) and 
are used to estimate the downlink Doppler correction factors. 


di = 


C<4 = 


/ (c - |v 8 / c -®(^ 

i/e )i)(c+iv/c-®(t: /e )D 


c 2 

ydss — Ofj-dss'j 

x d8S - 0 (f d88 ) _ x 8 / c -®(<- /c ) 

C 

|x ds8— ® (t dss ) - x 8 / c -®(fj /c )| 

V s/C — © (t S / C ) 

x d88 -®(tf sa ) - x s / c “®(t- /c ) 

C 

| x dss -0( f dss) _ x 8 / c -®(t’ /C )| 

j[c- |v d9S -°(tf 3 )|][c+ |v d88 -®(tf 3 )|] 


(16) 


(17) 


(18) 


(19) 


The gravitational redshift correction factors that account for the frequency shifts due to the masses of 
the Sun and the Earth, evaluated at t,, are given by 

1 1 

|x ds8 ~®| |x s / c_0 |. 

1 1 ' 

| x dss-©| | x® / c ©| 

The radiated 2.3-GHz frequency of the spacecraft is computed from the received frequency of the first 
data point as follows: 


j _ G M® 
d © “ c 2 


. _ G M 0 

d e ~ 


fs/c(t i) 


f'ih) 

(d\d2/d3di) + d© + d© 


+ f ; 


spin 
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where the Doppler and redshift terms are evaluated at t\. The estimated sky frequencies (evaluated at U) 
are computed from this transmitted spacecraft frequency at the first data point (evaluated at t{) using 
the above correction factors (evaluated at f ,) as follows: 


f S (ti) = fs/c(t l) 


did? 

d%d 4 


+ d© + d e 


For the one-way case, the residuals are defined such that the first data point is anchored to zero. Given 
the set of observable frequencies output from the OBSERVE F50 file, the estimated frequency for 

the first data point is, by definition, 


f s (h) = f a (h) 


The residual frequencies are computed by differencing the model from the observable as follows: 


A /'(fi) = /*(«<) -/'(ti) 


The residual frequencies and the spacecraft transmitted frequency are passed to STBLTY in the F52 file 
for further processing. 

3. Two-Way Data Processing. This section discusses the computations of two-way Doppler 
frequencies and residuals. The downlink Doppler corrections [Eqs. (16)-(19)] are evaluated at the time 
tag of the received signal frequency, tf 3S , as described in the previous section. In addition, the uplink 
vectors from the F49 file are interpolated to the time of uplink signal transmission, t“ p , and the Doppler 
correction terms for the uplink signal leg are similarly computed as follows: 


u \ (^z) — 'y 

/(c - |v dss_ 0 { f“ p ) | )( c + |v ds *-©(t" p ) | ) 

(20) 


c 2 

y-2{U) = 1 

v 8 / c -©(t*/ c ) 

x s / c - 0 (f’ /c ) - x d »»-0(t“ p ) 

(21) 

c 

| x »/c-O( t Vc ) _ x ds.-0(^P ) | 

U 3 (ti) = 1 

v d88 -0(t“ p ) 

x*/ c -®(t* /c ) - x d “-®(t“ p ) 

(22) 

c 

|x«/c-0( ( V c ) _ x d M -0 (t «P)| 

u 4 (ti) = y 

1 (c - |v 8 / c -©(fV' 

c )|)( C +|v»/ c -0(^ /c )|) 

(23) 


c 2 


The radiated frequency at the time of signal transmission, t“ p , is obtained from the uplink ramps in the 
PRF file, interpolated at t“ p , and then multiplied up to sky frequency, /*'(<“ p ). The estimated 2.3-GHz 
received frequency is computed as follows: 


/'(«<) 


rst/.up\ ^40 UjU2 _ , s /a 

u U ' 221 u 3 u 4 /jptn 


d\d2 
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The 2.3-GHz residual frequency is computed from the observed and estimated received frequencies at 
time ti as follows: 


A f s (ti) = f s (ti)-f s (ti) 


Likewise, for 8.4 GHz, 

did,2 

d 3 d^ 


f x (U) = 


rst (f u P\ 11 U\U2 _ ~ s /x 

Ju } 3 221 u 3 u 4 1 


spin 


a nu) = nu) - f x (u) 


4. Correction for Off-Axis Antenna on a Spinning Spacecraft. A sinusoidal signature due to 
an off-axis antenna on a spinning spacecraft can be removed by fitting a three-parameter model to the 
bias-corrected residuals: 


A / (ti) = C\ sin(C2ti + C 3 ) 


where ti is the time tag of the data point. The spacecraft transmitted frequency is corrected so that it 
is referenced to the phase center of the spacecraft spin. This is done by removing the value at the first 
data point: 


fs/c(ti) = fs/c(ti) - Cl sin(C 3 ) 


The residuals are adjusted accordingly: 


A f a c {t t ) = Af s (t l ) + Ci sin(C 2 fi + C 3 ) - C x sin(C 2 fi + C 3 ) 

The spacecraft transmitted frequency, f 3 / c , is passed to the program STBLTY in the F52 file along with 
other header information, observed frequencies, and residuals. 

E. STBLTY— Media Correction, Display Plots, and Output File 

STBLTY is the final program that applies additional corrections to the data, writes out records to a 
database file for USO and gravitational redshift analysis, and produces plots. At this point, the residual 
frequencies Af s (ti) and/or A / x (ti) for L data points i = 1,2, • • • , L are read in from the RESIDUAL F52 
output file. If both bands are available, STBLTY then combines the observable frequencies to produce 
the differential data type: 

A f sx (t l ) = f s (ti)-^f x {t i ) (24) 

For one-way data, the spacecraft transmitted frequency referenced to the first data point, f 9 / c (t 1 ), is 
passed in the header. 

1. Ionosphere/Charged-Particle Calibration. The charged particle effect can be “removed” 
from the radio science data in one of two ways: (1) When dual-frequency data are available, a differential 
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correction computed from the differential data type, Eq. (24), is applied to the individual frequency bands 
or (2) a TSAC-supplied polynomial describing the ionosphere path delay signature can be converted into 
frequency and removed from the frequency residuals. The TSAC correction removes ionospheric effects 
for either dual- or single-band downlink passes. 

For the case where two simultaneous downlink frequency channels exist, the charged particle can be 
removed using the differential frequency type f s - 3/711. given in Eq. (24). Given that the frequency 
residuals to be corrected are A f 9 (U) for 2.3 GHz and Af x (t t ) for 8.4 GHz, the corrected residuals are 


Af!(U) = A f(U) - — Af*{ti) 


for a one-way 2.3-GHz downlink; 


A/c(*0 = Af x (U) - f^Ar( tl ) 


for a one-way 8.4-GHz downlink; 


A fc(U) = A f'(U) - 


240 121 121 

221 112 + Il2 


A nu) 


for a two-way 2.3-GHz uplink/2. 3-GHz downlink; and 


A f!(ti) = A r(U) 


240 121 33 

22l 112 + U2 


Af sx (ti) 


for a two-way 2.3-GHz uplink/8. 4-GHz downlink, where the subscript c denotes the corrected value of 
the residual frequency. 

The TSAC group produces and delivers ionosphere calibration files to the RSSG, where they are 
archived and sometimes delivered to the radio science experimenters. The files contain polynomial coeffi- 
cients describing the mapped ionosphere path delay in meters (referenced to a 2295-MHz carrier) for each 
DSCC and for each time period that overlaps a tracking pass for a given spacecraft. TSAC determines 
these coefficients from observations of Earth-orbiting satellites at each complex and maps to the line of 
sight of the spacecraft. 


The TSAC coefficients and time tags are converted to frequency corrections in Hz (referenced to 
2.3 GHz) by STBLTY. Given the start time, (4), and end time, (£/), of the interval for which the 
polynomial is defined, and the polynomial coefficients (a,, i = 0, ■ • , 5), the ionosphere path delays to the 
spacecraft at regularly spaced intervals of time tj are computed as follows: 

5 

Api 0n {tj) = ^ ^ 

i=0 


where 


x 1 = 2 


tj 4 
tf ~ h 


- 1 
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The path delays computed at time offsets ±T/2 about tj are differenced and scaled to yield the one-way 
2.3-GHz ionosphere frequency correction at t 3 : 


/ a [Apion {tj + (T/2)) — A pion (tj (T*/2))] 

/ ion ) — C J> 

Figure 2 is an example of a typical path delay profile derived from a TSAC polynomial. Figure 3 is the 
corresponding ionosphere frequency profile at 2.3 GHz. 

The set of time tags, tj, and corrections, fi on {tj), that overlaps the tracking pass is used to perform 
the calibration only if dual-frequency charged-particle calibration was not previously performed. Given 
the frequency residual, A/ S (f z ) or A f x (U), and the above ionosphere correction interpolated to time 
ti, f s n (ti) , the corrected residuals are then computed as follows: For a one-way 2.3-GHz downlink, 


A r c (U) = A /'(«*) - f* on (t t ) 


for a one-way 8.4-GHz downlink, 


A/ C x (fi) = a nu) - ~f: on (u) 
for a two-way 2.3-GHz uplink/2. 3-GHz downlink, 

240 

A f' c (U) = A nu) - f° on (U) - — - tr) 

and for a two-way 2.3-GHz uplink/8.4-GHz downlink, 

3 11 240 

A /*(**) = A f X (ti) - -/^ n (ti) - y - tr) 


where t T is the round-trip light time. 

For the case of differential / 2 3 - 3/ x /ll two-way data, the contribution of the common 2.3-GHz uplink 
will be zero after differencing the two downlink channels. The remaining charged-particle contribution 
in the A f 3X data type (SX) will be due to the downlink only. Therefore, this data type is corrected as 

follows: 


A /“(ti) = A r x (U) 


—f s 
1 \ 1% ° 




(25) 


2. Troposphere Calibration. The troposphere calibration is performed after the ionosphere cali- 
bration so that the remaining time-dependent signature of the residuals over elevation angle more closely 
follows the troposphere ray-path function. The troposphere can be removed using either a model with a 
specified zenith troposphere path delay (default 2.1-m) or the zenith path delay can be fit and removed 

from the data if there is sufficient arc over elevation angle. The ray-path mapping function used is the 
CfA (Center for Astrophysics) mapping function [14]: 
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1 


S \ne(t t ) + (6i/(tane(i) + (f>2/(sine(ii) + 63))) 
where e(t») is the elevation angle and the coefficients b k are functions of meteorological parameters, given 

by 


= 0.0002723[1 + (2.642 x 1CT%) - (6-400 x 10“%) + (1-337 x 10" 2 T o ) - (8.550 x 10 2 a) 

- (2.456 x 10 ~ 2 h 2 )} 

b 2 = 0.0004703(1 + (2.832 x 10" 5 po) + (6.799 x 10-%) + (7.563 x 10"%) - (7.390 x 10- 2 a) 

- (2.961 x 1( r 2 h 2 )} 

63 = -0.0090 

Reasonable values of surface pressure (po — 900 mbar), partial water vapor pressure (eo 5 mbar), 
surface temperature (T 0 = 292 K), tropopause altitude (h 2 = 12.2 km), and dry model parameter 
(a = 5.0) are used in the computation of these coefficients. This model is in agreement with predictions 
made by Lanyi [15] to the 1-cm level down to a 5-deg elevation angle. 

The one-way troposphere correction to the observed 2.3-GHz frequency residual is given by 


flrJto) = ~ 


f s {tj)Ap z 

cT 


where f s {U) is the downlink sky frequency, A p z is the zenith troposphere path delay, c is the velocity 
of light, and T/2 is a chosen time offset about t % . For the cases of two-way 2.3-GHz uplink/2.3-GHz 
downlink and 2.3-GHz uplink/8. 4-GHz downlink, the corrections are 


ftropifi) — 


ftrop(ti ) ~ 


-/%,)Apz 

cT 


cT 



where t r is the round-trip light time. The one-way frequency residuals are adjusted accordingly: 


A r c (U) = A f‘(U) - f t s rop (U) + f t ‘ rop (t 1) 


A/*(t.) = A r(U) - f t % op {U) + f? rop (t 1) 


effectively removing the signature of the troposphere relative to the first data point, because fs/c{t 1) 
was evaluated here. A bias correction evaluated from the average of the corrected residuals is applied 
to the transmitted spacecraft frequency to refer it to the center of mass of the residuals. The spacecraft 
transmitted frequency at 2.3 GHz is corrected as follows: 


/c%/c(W) = f! /c (tl) + ft 3 rop(tl) + A n 


where the last term is the average of the corrected residuals over the data span and t mi d is the time tag 
of the midpoint of the data acquisition interval. 
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The two-way frequency residuals are corrected for troposphere as follows: 


A/‘(t t ) = A f(U) ~ ftropiU) 

A f X c{U) = A /«(«<) - !? Top (U) 


III. STBLTY Performance Analysis Results and Discussion 

The following sections will discuss examples using STBLTY on Galileo and Ulysses flight data. 

A. Galileo Results 

Examples of running STBLTY on Galileo radio science data will be presented in the following two 
sections. 

1. One-Way Data. The results of routine processing of Galileo USO data using STBLTY will be 
briefly discussed in this section. For further details on this analysis of closed-loop one-way USO data 
acquired during the first 2 years of cruise (1989-1991), the reader is referred to [3]. 

Due to the undeployed high-gain antenna (HGA), only 2.3-GHz data from the LGA antennas were 
available. The passes were conducted weekly on the average, were typically 2 hours in duration, and were 
centered about meridian crossing where the troposphere effect was minimal. For troposphere calibration, 
it was sufficient to remove a model (versus fit and remove) using a default zenith path delay. The error 
incurred by not performing an ionosphere calibration is expected to be at the mHz level or below for a 
2-hour pass centered about meridian crossing at 2.3 GHz. 

The Allan deviation results were consistent with the values expected at the appropriate time intervals, 
given the known signal levels and noise mechanisms. At the short time intervals, the noise in the signal 
was dominated by thermal noise due to the low signal levels that were due to transmitting through 
the LGA. At the long time intervals, the noise in the signal was dominated by the USO itself. The 
transmitted USO-referenced spacecraft frequency followed a behavior that was consistent with known 
aging mechanisms of the USO, having an rms scatter of 17 mHz about a fitted aging model. 

2. Two-Way Data. The testing of different calibration schemes using STBLTY on Galileo two-way 
2.3-GHz data will be discussed in this section. STBLTY was run on three passes: (1) a DSS-63 tracking 
pass on day 91-123, (2) a DSS-63 pass on day 93-081, and (3) a DSS-14 pass on day 93-082. The latter 
two passes are from the 1993 Gravitational Wave Experiment. 

Table 2 displays the processing results for these three passes using different calibration schemes. The 
first column denotes the year and day number of the pass; the next column is the DSN station identifica- 
tion; the next column denotes the ionosphere calibration employed (see the key for an explanation of the 
codes); the next column denotes the troposphere calibration employed (see the key for an explanation of 
the codes); the next column is the zenith path delay and uncertainty in the troposphere calibration (if 
the code in the TRP column is ll F,” this value is from the least-square fit); the next column is the slope 
from a linear fit of the residuals after calibrations have been performed; and the last four columns are 
the Allan deviations of the postcalibration residuals at 1, 10, 100, and 1000 s. 

Four different calibration schemes were employed for each pass. The first run (ION:T, TRP:F) involved 
calibrating the ionosphere utilizing TSAC polynomial coefficients and calibrating the troposphere by 
fitting the model over the elevation angle for a zenith path delay and bias. The second run (ION:N, 
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TRP:F) involved no ionosphere calibration and calibrating the troposphere as described in the first 
scheme. The third run (ION:T, TRP:M) involved calibrating the ionosphere using TSAC polynomial 
coefficients and calibrating the troposphere by removing the model using the default zenith path delay of 
2.1 m (no fitting performed). The last run (ION:N, TRP:N) did not employ any media calibration. 

The zenith path delay for the TSAC ionosphere/troposphere fit run (T,F) for passes 91-123 and 93-081 
are in reasonable agreement with the expected 2.1-m zenith path delay. The zenith path delay fit from the 
data when no ionosphere is removed (N,F) is consistently lower than those of the TSAC calibrated runs 
(T,F) and (T,M). This suggests that the TSAC calibration does contribute to removing nontropospheric 
trends from the data. 

The linear slope of the residuals tends to increase in magnitude as fewer calibrations are performed. 
The magnitude of the linear slope tends to be smaller for the runs involving both troposphere fit and 
TSAC ionosphere calibration (T,F) and larger for the run involving no calibrations (N,N). 

The Allan deviations at the time intervals of 1, 10 and 100 s do not differ significantly between the 
different schemes over each given pass. The Allan deviation at 1000 s for the three runs involving a tropo- 
sphere calibration (T,F), (N,F), and (T,M) are in reasonable agreement. The 1000-s Allan deviations are 
consistent with values predicted due to solar plasma at the known solar elongation angles using the model 
of Armstrong, Woo, and Estabrook [16]. However, the Allan deviations at 1000 s are significantly larger 
for the no-media calibration case (N,N), clearly showing that the removal of the systematic troposphere 
signature with elevation angle is important for obtaining meaningful results at this time interval. 

For Galileo at 2.3 GHz, trajectory errors of 0.018-mm/s in velocity translate to Allan deviations at 
1000 s of about 6 x 10~ 14 . This lies below the 1(T 13 level expected due to plasma (for two-way data) and 
the instability of the USO (for one-way data). 


Table 2. Galileo media calibration test results. 


Pass, 

yr-day 

DSS 

ION 

TRP 

Zenith path 
delay, 
m 

Linear 
slope, 
10 -7 Hz/s 

^v(l), 

ltr 14 

10“ 14 

Ml 00 ). 

IQ - 14 

CTy(lOOO), 

io -14 

91-123 

63 

T 

F 

1.96 ± 0.03 

0.12 

± 

0.29 

- 

357 

27.7 

24.8 

91-123 

63 

N 

F 

0.98 ± 0.03 

0.16 

± 

0.29 

— 

357 

27.6 

27.0 

91-123 

63 

T 

M 

2.1 

0.84 

± 

0.29 

— 

357 

27.7 

25.0 

91-123 

63 

N 

N 

— 

-5.01 

± 

0.29 

— 

357 

27.7 

38.2 

93-081 

63 

T 

F 

2.17 ± 0.03 

0.77 

± 

0.27 

1947 

206 

25.4 

17.9 

93-081 

63 

N 

F 

1.70 ± 0.04 

1.55 

± 

0.27 

1947 

206 

25.5 

23.3 

93-081 

63 

T 

M 

2.1 

0.58 

± 

0.27 

1947 

206 

25.4 

18.4 

93-081 

63 

N 

N 

— 

-3.55 

± 

0.27 

1947 

206 

26.1 

50.5 

93-082 

14 

T 

F 

1.11 ± 0.04 

0.39 

± 

0.32 

1808 

194 

26.9 

14.5 

93-082 

14 

N 

F 

0.80 ± 0.05 

0.75 

± 

0.31 

1808 

195 

26.9 

17.9 

93-082 

14 

T 

M 

2.1 

2.73 

± 

0.32 

1808 

195 

26.9 

13.4 

93-082 

14 

N 

N 

— 

-1.14 

± 

0.32 

1808 

195 

27.1 

29.1 


T = TSAC ionosphere calibration. 

D = Dual-frequency ( f s - 3/ x /ll) charged-particle calibration. 

N = No calibration performed. 

F = Troposphere calibration (model removed using fit zenith path delay). 

M = Troposphere calibration (model removed using default zenith path delay). 
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B. Ulysses Results 

For Ulysses, two sets of two-way data were analyzed: (1) closed-loop data from a dual-band pass 
from the first opposition 91-004 (DSS 14) and (2) open-loop data from a dual-band pass from the second 
opposition 92-064 (DSS 14). Table 3 displays the results for these data sets and for each of the applicable 
data types within each pass. For the single-band data (2.3 and 8.4 GHz), the following runs were 
performed: (1) TSAC ionosphere and troposphere fit (T,F), (2) dual-frequency calibration of charged 
particles, and troposphere fit (D,F), (3) no ionosphere/charged-particle and troposphere fit (N,F), and 
(4) TSAC ionosphere and troposphere model removed using the default zenith path delay (T,M). The 
no-calibration case (N,N) is not presented in this table since the (N,N) Galileo runs (Section III. A) clearly 
illustrate the degradation when not applying a media correction. The description of Table 3 is the same 
as that of Table 2 in the proceeding section except that the data type code (S, X, or SX) follows the 
year-day number in the first column. 


Figure 4 is a plot of the 8.4-GHz residuals for pass 91-004 for the fully calibrated (T,F) run. Figure 5 
is a plot of the averaged residuals where every 60 points have been averaged to allow long period trends 
to be easily examined. Figure 6 is a plot of the reconstructed phase for this pass, and Fig. 7 is a plot of 
the corresponding Allan deviation. Figure 8 is a plot of the averaged residuals prior to the removal of the 
troposphere, thus illustrating the signature of an unmodeled troposphere. Figure 9 displays the elevation 
profile over the pass and the troposphere model that was fit from the data of Fig. 8 and then removed 
from the residuals, resulting in the plots displayed in Figs. 4-7. Figure 10 displays the TSAC ionosphere 
correction at 2.3 GHz for this pass. 


Table 3. Ulysses media calibration test results. 


Pass, 

yr-day 

DSS 

ION 

TRP 

Zenith path 
delay, 
m 

Linear 
slope, 
10~ 7 Hz/s 

*#(1). 
10“ 14 

*1,(10), 

10-14 

*#(100), 
10— 14 

(Ty(1000), 

io- 14 

91-004 S 

14 

T 

F 

1.50 ± 0.06 

-0.001 ± 0.71 

1080 

116 

15.4 

4.8 

91-004 S 

14 

D 

F 

1.90 ± 0.08 

0.08 ± 1.71 

2624 

276 

32.5 

2.4 

91-004 S 

14 

N 

F 

0.44 ± 0.06 

0.09 ± 0.71 

1081 

116 

15.4 

6.2 

91-004 S 

14 

T 

M 

2.1 

1.35 ± 0.71 

1081 

116 

15.4 

6.7 

91-004 X 

14 

T 

F 

1.69 ± 0.03 

0.29 ± 1.65 

621 

66 

9.3 

1.9 

91-004 X 

14 

D 

F 

1.40 ± 0.04 

0.65 ± 1.83 

756 

80 

12.8 

3.7 

91-004 X 

14 

N 

F 

1.05 ± 0.04 

0.42 ± 1.65 

621 

66 

9.3 

2.8 

91-004 X 

14 

T 

M 

2.1 

3.51 ± 1.65 

621 

66 

9.3 

3.6 

91-004 SX 

14 

T 

N 

— 

0.66 ± 0.67 

1107 

116 

13.3 

3.3 

91-004 SX 

14 

N 

N 

— 

1.51 ± 0.67 

1107 

116 

13.4 

5.2 

92-064 S 

14 

T 

F 

1.54 ± 0.05 

0.19 ± 0.12 

335 

46 

16.7 

6.7 

92-064 S 

14 

D 

F 

1.65 ± 0.04 

-0.004 ± 0.12 

384 

45 

16.4 

7.4 

92-064 S 

14 

N 

F 

1.17 ± 0.04 

0.13 ± 0.12 

335 

46 

16.7 

5.9 

92-064 S 

14 

T 

M 

2.1 

1.26 ± 0.12 

335 

46 

16.7 

6.6 

92-064 X 

14 

T 

F 

1.43 ± 0.03 

0.12 ± 0.25 

147 

30 

13.6 

5.6 

92-064 X 

14 

D 

F 

1.46 ± 0.04 

0.15 ± 0.27 

187 

31 

13.8 

6.1 

92-064 X 

14 

N 

F 

1.38 ± 0.04 

0.23 ± 0.25 

147 

30 

13.6 

5.7 

92-064 X 

14 

T 

M 

2.1 

4.63 ± 0.24 

147 

30 

13.6 

6.2 

92-064 SX 

14 

T 

N 

— 

0.15 ± 0.09 

303 

32 

9.4 

4.1 

92-064 SX 

14 

N 

N 

— 

0.49 ± 0.09 

303 

32 

9.4 

4.1 
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Fig. 5. An example of 8.4-GHz frequency residuals from Fig. 4 averaged every 60 s from 
a Ulysses pass conducted at DSS 14 on January 4, 1991. 
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Fig. 10. Ionosphere frequency correction at 2.3 GHz, derived from a TSAC polynomial 
for the DSs-14 Ulysses pass of January 4, 1991. 


The fitted zenith path delays for case (T,F) appear reasonable but are somewhat low compared to 
the expected 2.1-m values. Unmodeled trends, which blend with the signature of the troposphere, are 
assumed to exist in the data and thus influence the fit. 

The linear slope from a straight line least-squares fit of the postcalibrated residuals tends to be of 
insignificant or of marginal magnitude when a troposphere fit is performed. When a model troposphere 
is removed, the fitted linear slopes of the residuals are significant. This suggests that significant unknown 
trends may be present in the data set that are absorbed by the troposphere fit or that the troposphere 
model using the default zenith path delay is not sufficient. 

The Allan deviations at 1, 10 and 100 s do not change significantly between the different runs for a 
given band, except for the dual- frequency charged-particle calibrated runs (D in the ION column), where 
significant system noise at these time intervals is introduced into the data. At 1000 s, dual-frequency 
charged-particle calibration (D,F) results in the lowest Allan deviations at 2.3 GHz for the 91-004 pass 
and no significant changes at 2.3 GHz for pass 92-064. At 8.4 GHz, the TSAC calibrated run for pass 
91-004 produced the lowest Allan deviation, 1.9 x 10 For the 92-064 DSS-14 pass, no significant 
improvement was noted in the 8.4-GHz 1000-s Allan deviation between the different runs. The 91-004 
Allan deviations agree with independently measured values [4] . 

For the f s - 3/ x /ll data type (SX), the following runs were performed: (1) TSAC correction 
[Eq. (25)] (T) and no troposphere (N) and (2) no ionosphere calibration (N) and no troposphere cali- 
bration (N). Since the difference frequency eliminates all nondispersive effects, troposphere calibration 
is not performed. For 91-004 SX and 92-064 SX, the slopes of the calibrated runs were of insignificant 
or marginal magnitudes, suggesting the effectiveness of using TSAC polynomials to remove long period 
trends from the SX data type and illustrating the dominance of the time variability of the ionosphere 
over that of the solar plasma during these observations near solar opposition. The Allan deviations 
do not change significantly between runs for the SX data type, except that the 1000-s value for pass 
91-004 SX shows an improvement from 5.2 x 10“ 14 to 3.3 x 10~ 14 when applying the TSAC calibration. 
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Figure 11 displays the averaged residuals for this pass without applying TSAC (N,N). Here, all nondis- 
persive effects should cancel out, and the resulting signature is presumably due to dispersive sources. The 
corresponding TSAC S-X frequency calibration for this pass is shown in Fig. 12. The averaged residuals 
after applying the TSAC calibration (T,N) are shown in Fig. 13. The resulting signature is considerably 
flatter than that of the non-TSAC calibration run (N,N) of Fig. 11, but a long period trend of order 
0.5-mHz over this pass is apparent upon closer inspection. This is consistent with the observation of 
Bertotti et al. [4] that the error of the TSAC calibration is not negligible and can be as large as 0.5 mHz. 
The long period trends on the order of 0.5 mHz that remain in the data after applying media calibrations 
are consistent with known uncertainties in the calibration schemes. The zenith path delays for the TSAC 
ionosphere can be in error by as much as 75 cm for day passes and 15 cm for night passes. The troposphere 
can be in error by 5 cm. 9 * An overall 15-cm zenith path delay error (for these nighttime solar opposition 
passes) can map into variations at the 1-mHz level over a pass, consistent with what is observed in these 
data. 

For Ulysses at 8.4 GHz, trajectory errors of 0.012 mm/s in velocity at 1000 s translate to Allan 
deviations at about 4 x 10“ 14 , which is comparable with the observed values. These are also comparable 
with expected noise due to media. After removal of the ionosphere and troposphere calibration models, 
significant random fluctuations in the data remain. These fluctuations have different characteristics with 
different frequencies, solar elongations, local weather conditions, and local ionospheric variability. These 
are usually masked by system noise at the lower time intervals, whose magnitudes depend upon received 
signal levels and received noise power of a given spacecraft. Expected levels of these fluctuations at 1000 s 
for Ulysses at 8.4 GHz are about 3 x 10~ 14 for the troposphere, an upper bound of 3 x 10" 14 for plasma 
noise, and 2.5 x 10“ 14 for the system noise of the closed-loop receiver [4]. 

Figure 14 displays the 8.4-GHz Allan deviation curves for a set of different calibration schemes for the 
91-004 DSS-14 pass. For the case of two-way data over a sufficient elevation angle arc, the best estimate 
of system stability is obtained by applying a TSAC ionosphere correction and by fitting and removing a 
troposphere model. 


IV. Conclusion 

A description of the STBLTY program set has been presented with details on model implementation 
and residual analysis. Examples of the use of this program set on Galileo and Ulysses data were also 
presented. The frequency of a spacecraft signal received on the ground is the observable data type. All 
known effects are removed from this observable to produce residuals that can be inspected to evaluate 
performance. Among the errors in the residuals are those due to white noise, media, and trajectory. For 
a sufficiently high SNR, random media effects appear to be the limiting error source at the 1000-s time 
scale for both Galileo and Ulysses coherent data sets. 

Future upgrades to STBLTY will focus on modeling the gravity fields of planets in order to remove 
these effects from the data during flybys or orbital operations. The first use for this capability will involve 
the Galileo Jupiter orbital operations when Galileo is expected to tour the Jovian system starting in late 
1995. Other future upgrades include modifying the code to process 32.0-GHz downlink data such as that 
expected for the Cassini radio science experiments, to process coherent downlink data that involve 8.4- 
or 32.0-GHz (Ka-band) uplink signals, converting the code to manipulate vectors in the J2000 system, 
and incorporating a planetary atmospheric model for occultation data analysis. The program set was 
recently modified to process three-way data (implemented in a developmental version). 


9 T. McElrath, personal communication, Navigation Systems Section, Jet Propulsion Laboratory, Pasadena, California, 

1994. 
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FREQUENCY, Hz AVERAGED FREQUENCY RESIDUALS, Hz 



Fig. 11. Averaged differenced frequency residuals (dual-band) for the DSS-14 
Ulysses pass of January 4, 1991. 



Fig. 12. Ionosphere frequency correction for dual-band derived from a TSAC 
polynomial for the DSS-14 Ulysses pass of January 4, 1991. 
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AVERAGED FREQUENCY RESIDUALS, Hz 




Fig. 14. Allan deviation curves for a set of different calibration schemes for 
the DSS-1 4 Ulysses pass of January 4, 1 991 . 
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Appendix 

Formulation of Precession and Nutation Matrices 

The following definitions for system rotation matrices and their derivatives apply when defining the 
precession and nutation rotation matrices and their derivatives: 

[i o o 1 d r 0 o o] 

{9)x= 0 cos 9 sin 9 — (#)* = 0 — sin (9 cos# &=(--$) 9 

.0 -sin# cos# J [0 -cos# -sin# ' 2 x 

cos# sin# Ol , f-sin# cos# 0] 

(9) z = -sin# cos# 0 — (9) z = -cos# -sin# 0 6 = (- - 6) 9 

.0 0 1J dt 0 0 Oj ' 2 z 

I. Precession Rotation Formulation 

The precession matrix, P, is the product of the following three system rotations using the definitions 
given above: 

^ = (2 ~ a 1950^ (^— ~ ^1950 j (”Ai95o) 2 

where the precession rotation angles in terms of Julian centuries since 1950, T, are given by 

^1950 = [89.9999986565 - 0.64027801T - (3.042075 x 10 _4 )T 2 - (5 0837 x 10~ 6 )T 3 1-^- 

v } j 180 

^1950 = [89.9999988317 - 0.5567500297T 4- (1.185607 x 10“ 4 )r 2 + (1 16119 x lO" 5 ^ 3 ] — 

; J 180 

a i960 = [(-1.3435 x 10~ 6 ) - 0.6402780091T - (8.39481 x 10~ 5 )T 2 - (0.50003 x 10 _5 )T 3 1 — 

; J 180 

The derivative of the precession matrix is as follows: 

P = (7T - a 195 o)^ - £ 195 o) (-£ 1950 ), “^ a 1950 

+ (2 ” Q 195oj^ (tt - (5i95o) x (-Ai95o) 2 — <$1950 

+ (I ■ ai95 °), (I - 6 ' 950 ) x (f - Al950 ) 2 f^r Al950 


where 
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^Aigso = [-0.64027801 - (2 x 3.042075) x 10" 4 T - (3 x 5.0837) x 10-®T 2 ]^ 

^1950 = [-0.5567500297 + (2 x 1.185607) x 10 _4 T + (3 x 1.16119) x 10 ~ 5;r2 ]^ 

^aigso = [-0.6402780091 - (2 x 8.39481) x 10~ 5 T - (3 x 0.50003) x lO^T 2 ]^ 

II. Nutation Rotation Formulation 

The nutation rotation matrix is the product of the following three rotations using the matrix rotation 
definitions given above: 

N = (-e m ) x (A^(e m + Ae) x 

where the nutation angles are the IAU 1980 model values defined in the text. The derivative of the 
nutation rotation matrix is thus given as 

N = -e(e m + |) (A¥),(e m + Ae) I + A*(-e m )x(-Atf + ^) i 
x (c m + Ae) x + e(-e m )x(A^) z (-e m - Ae + 
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Analysis of a Microstrip Reflectarray Antenna 
for Microspacecraft Applications 

J. Huang 

Spacecraft Telecommunications Equipment Section 


A microstrip reflectarray is a Hat reflector antenna that can be mounted confor- 
mally onto a spacecraft’s outside structure without consuming a significant amount 
of spacecraft volume and mass. For large apertures (2 m or larger), the antenna’s 
reflecting surface, being flat, can be more easily and reliably deployed than a curved 
parabolic reflector. This article presents the study results on a microstrip reflect- 
array with circular polarization. Its efficiency and bandwidth characteristics are 
analyzed. Numerous advantages of this antenna system are discussed. Three new 
concepts using this microstrip reflectarray are also proposed. 


I. Introduction 

JPL is currently developing microspacecraft technologies for future deep space missions in order to 
meet NASA’s goal of having small, efficient, and inexpensive spacecraft. The microspacecraft, having sizes 
on the order of one-half meter, will certainly require components that are small both in size and mass. 
High-gain antennas are one part of the telecommunications equipment that warrants attention since they 
generally require a significant amount of real estate and mass. The conventional high-gain antennas most 
often used are parabolic reflectors. Although they are efficient radiators, parabolic reflectors are generally 
bulky in size and large in mass, due to their curved reflecting surfaces. As a result, a flat reflector called 
a microstrip reflectarray [1,2] is being proposed as a future candidate high-gain antenna. It is well known 
that when a required antenna gain is given at a particular frequency, the antenna aperture size is more or 
less fixed. The only significant size reduction that may be achieved for an antenna is its profile thickness. 
The flat-plate microstrip reflectarray offers such an advantage of profile size reduction as compared to a 
conventional parabolic reflector. 

The reflectarray antenna as shown in Fig. 1 represents World War II technology [3], However, the 
low-profile printed reflectarray is a fairly new concept. The microstrip reflectarray, being one of the 
printed low-profile antenna technologies, consists of a very thin, flat reflecting surface and an illuminating 
feed, as shown in Fig. 2. On the reflecting surface, there are many isolated microstrip patch elements 
without any power division network. To each patch element is attached a short segment of phase-adjusting 
transmission line to compensate for the phase delay over the path from the illuminating feed. Because of 
the phase adjustment capability of the patch elements, the reflecting surface can be flat or conformal to 
its mounting structure and still maintain a constant phase aperture field. A detailed description of this 
antenna concept, as well as its advantages, is given in the following sections. Theoretical analysis of the 
antenna performance parameters, such as radiation pattern, efficiency, and bandwidth, is also presented. 
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Fig. 1. Reflectarray configuration: (a) three-dimensional view and (b) two-dimensional view. 
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Fig 2. Flat-plate microstrip reflectarray with identical patches but different-length 
microstrip transmission lines. 


II. Description of the Microstrip Reflectarray 

When many antenna elements with open- or short-circuited terminations are arranged in a planar 
aperture and are illuminated by a feed antenna, as shown in Fig. 1, these elements will reradiate their 
illuminated energy into space. The total reradiated energy will be noncophasal if all the elements and 
their terminations are identical. This is because the fields that propagate to the elements from the 
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feed have different path lengths, Si, S 2 , ■ • ■ , S^, as shown in Fig. 1(b), and thus form different phases. 
However, if each element’s phase is adjusted to compensate for these different path lengths, the total 
reradiated field can be made cophasal and concentrated toward a specific direction. The array antenna 
formed by the above concept is named the reflectarray and was introduced [3] many decades ago using a 
horn, dipole, open-ended waveguide, etc., for the element. Since these elements are large in size at lower 
microwave frequencies and many elements are needed in a reflectarray in order for it to be efficient, the 
earlier reflectarray antennas were bulky in size and heavy in weight. Due to the recent advancement of 
the lightweight and low-profile printed antennas, such as the microstrip patch, the printed reflectarray 
becomes physically more realizable and attractive. Several different versions of the printed reflectarray 
have recently been developed. One version, shown in Fig. 2, has all identical microstrip patch elements 
but with different-length microstrip transmission lines. This microstrip reflectarray was first patented by 
Munson and Haddad [4] and then openly published and analyzed by Huang [2] and Metzler [5]. Litva, 
Zhuang, and Wu [6], as well as Chang and Huang [7], have demonstrated the concept via hardware 
development. Malibu Research Center developed a different reflectarray [8] with printed dipoles having 
different sizes and no phase delay lines. The required path delay-compensating phases from the elements 
are achieved primarily via the differing lengths of the dipoles. Different lengths yield different input 
impedances (complex quantity) at a particular frequency, which in turn give different phases. Targonski 
and Pozar [9] developed a microstrip reflectarray with patches having different sizes and no delay lines. 
It is expected that the reflectarray with printed dipoles or patches having different sizes will not radiate 
as efficiently as the identical-patch reflectarray with different-length delay lines. This is because, for a 
particular frequency, there is only one optimal size of the resonant structure (dipole or patch) to reflect or 
reradiate energy; other sizes will result in lower amplitudes. This is similar to the phenomenon observed 
with frequency selective surfaces (FSS), where only a narrow spectrum will yield total reflection. 

By reason of its expected higher efficiency, the microstrip reflectarray with patches having identical 
sizes and different phase-delay lines is proposed and studied here. This flat reflector antenna, as shown 
in Fig. 2, is composed of a thin (<0.03 wavelength) slab of dielectric material having one side completely 
covered with a thin layer of metal (serving as a ground plane) and the other side etched with many 
identical metallic microstrip patches. A feed antenna, located at an optimally designed distance from 
the array elements, will effectively illuminate all the patches. The size of each patch, which can be 
rectangular, square, or circular, is designed to resonate at the frequency of the feed antenna. A short 
transmission line is connected to each patch at one end with the other end of the line either open or 
short circuited. To generate circular polarization, two equal-length transmission lines are needed to be 
orthogonally connected to each square or circular patch, as shown in Fig. 3. In this case, the feed has to 
be circularly polarized. The short transmission line connected to each patch can be either a microstrip 
line etched on the same side of the patch or a stripline sandwiched in an additional layered structure 
placed behind the patch’s ground plane. The advantage of the microstrip line is ease of fabrication, 
while that of the stripline is minimum interference to the patch’s radiation. When the radiation field 
of the feed antenna (in transmit mode) strikes each patch, the received resonant field of the patch will 
travel through its connected transmission line and be reflected by its open- or short-circuited termination 
and then reradiated by the patch into space. Thus, all the microstrip patches behave as reradiators, 
while the short transmission lines serve as phase delay lines. The lengths of these transmission lines are 
intentionally made different for differently located patches to compensate for the path delay differences 
from the feed antenna. With proper design and calibration of these line lengths, the reradiated fields from 
all the patches can be made cophasal in the broadside direction. Also, by redesigning the line lengths, 
the main beam can be directed toward other directions at large angles (up to more than 50 deg) from 
the broadside direction. Since the required phase changes for all the elements are between 0 and 360 deg, 
the maximum length needed for the transmission line is only one-half wavelength. Consequently, the 
insertion loss associated with these short lines will be insignificantly small. The transmission line should 
be impedance matched to the patch radiator, which can be done using a quarter-wave-long impedance 
matching section. Because this is a phase-delay approach and not a time delay, as frequency changes, 
a phase excursive error will occur, especially for the outer elements of the array (assuming the feed is 
located at the center axis of the array). In other words, the phase will accumulate more error for the 
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outer elements as the frequency changes, and this limits the bandwidth performance of the reflectarray. 
This accumulated phase error can be reduced by using longer transmission lines for the center elements 
(time delay approach) or by using a larger f/D ratio, where / is the distance between the feed and the 
array center and D is the array diameter. 

Since the microstrip reflectarray does not require any power divider, its efficiency in a large array 
system is much higher than a conventional array having the same aperture size. One possible drawback 
is that, in addition to the reradiated fields from the patches, there will also be scattered field from the 
patches, reflected field from the ground plane (especially away from the resonant frequency of the patch), 
scattered field from the phase delay lines, and diffracted field from the edge of the reflectarray. These 
backscattered fields may increase the sidelobe level and possibly distort the main beam shape. However, 
because most of these scattered fields are noncophasal, and as long as the aperture directivity of the 
reflectarray is sufficiently higher (20 dB or more) than the feed directivity, the backscattered energy is 
generally small relative to the desired main beam. In other words, the microstrip reflectarray can be an 
efficient antenna system only if it has a large number of array elements (500 or more). 

(a) 


y 


£ 



Fig. 3. Reflectarray (a) square and (b) circular microstrip elements, for 
circular polarization. 


III. Radiation Analysis of the Microstrip Reflectarray 

When a rectangular planar array with M x N microstrip patch elements is nonuniformly illuminated 
by a low-gain feed at r/, as shown in Fig. 4, the reradiated field from the patches in an arbitrary direction, 
u, will be of the form 


M N 


E(u) 


EIM V rnn ’ f* / ^ ^ mn A t^o) ^Cp ^ j ko T' mn 


+ ja n 


= 1 n=l 


(1) 


where F is the feed pattern function, A is the pattern function of the microstrip patch, r mn is the 
position vector of the ranth patch, u 0 is the desired main-beam pointing direction, and o m n is the 
required transmission-line phase delay of the mnth element. The condition for the aperture distribution 
to be cophasal in the desired direction u Q is 

Omn |^| ^ mn ^ / [ T T mn ’ &o j = 2717T, 71 0, 1 , 2, ■ * * (2) 

For a circular aperture, which is desirable for better aperture efficiency as compared with a rectangular 
aperture, the summation signs in Eq. (1) can be truncated to 0 (no calculation) for patches located outside 
of the circular aperture. 
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Fig. 4. Coordinate system for reflectarray pattern analysis. 


Equation (1) gives a nonpolarized (scalar) field. For circularly polarized radiation, the field of Eq 
can be separated into 0 and 0 components: 


Eq = 

= E(u) 


)pi 

■0+(j 

2 ‘ P2 J 

P 2 


E$ - 

= E( u) 

[(/,'P,; 

I pi 

■0 + (j 

f - N 

2 • Vi t 

|P2 

■*] 


•( 1 ) 


(3a) 

(3b) 


where f\ and f 2 are the two orthogonal polarization vectors of the feed horn at the patch location, pi and 
p 2 are the two orthogonal patch polarization vectors, and 9 and <j> are the orthogonal far-field spherical 
coordinate unit vectors. For a right-hand circularly polarized feed, the far-field copolarized radiation from 
the reflectarray is 


Eco-pol — — f= {Eq + jE#) 


(4a) 


and the cross-polarized radiation is 


E x -pol - {Eq — jE 0) 


(4b) 


In Eq. (1), the feed pattern function F is modeled by a cos" $ function. For the pattern function A of 
the single microstrip patch, a simple closed-form model using the dual-slot theory [1] is employed. This 
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simple model, which is accurate enough for large array prediction, allows the computation time of many 
thousands of patch elements to become more realistic as compared with other more rigorous techniques. 
Consequently, mutual coupling effects between patches are not included. For a substrate thickness of less 
than 0.03 free-space wavelength, beam scan angles of less than 45 deg, and no extremely low sidelobe 
requirement, the mutual coupling effects can generally be neglected for the microstrip array. Experimental 
reports [10] demonstrated that, due to the low-profile nature of the microstrip antenna, mutual coupling 
has not been found to be a serious problem in most of the array applications. 

In addition to the reradiated field given in Eq. (1), there is a certain amount of backscattered field, 
such as the scattered field from the phase delay transmission lines, the patches, the ground plane, and the 
ground plane edges. However, as long as the reflectarray is designed with proper element spacing (avoiding 
the grating lobe condition) and the patch element is designed to resonate at the correct frequency, the 
backscattered fields can be minimized. Furthermore, if the directivity of the feed antenna is significantly 
lower (by at least 20 dB) than the directivity of the reflectarray aperture, the effect of the noncophasal 
backscattered field will generally be insignificant [1] when compared with the cophasal reradiated fields 
from the patches. To summarize, with proper design and large enough array size, Eq. (1) is an excellent 
approximation with which to efficiently calculate the far-held radiation pattern of the reflectarray. It 
should be able to predict the main beam and the first few sidelobes with good accuracy. 


IV. Antenna Efficiency Analysis 

The two primary factors that govern the efficiency of the microstrip reflectarray are very similar to 
those for the parabolic reflector. These are the aperture illumination efficiency and the feed spillover 
efficiency [11]. Other minor factors that contribute to the reflectarray efficiency are the patch element 
loss and the back scattered held loss. The aperture illumination efficiency is caused by the unequal 
illumination of the array aperture due to the feed’s tapered pattern. A uniformly illuminated aperture is 
dehned as having 100-percent illumination efficiency. The spillover efficiency is the ratio of the amount of 
feed energy that illuminates the entire array to the total amount of energy that is radiated by the feed, 
which includes the amount that spills outside the array aperture. It is clear that the illumination and 
the spillover efficiencies are complementary to each other. In other words, if one increases, the other will 
decrease, and vice versa. 

Let us define the total aperture efficiency (rj) as the product of the illumination (t)m) and spillover 
(tj s ) efficiencies: 


V = Vill ’ Vs 


( 5 ) 


From Silver [12] and Fig. 5, 


\jLCC -p^s\ 2 i/i 2 

*£ oCoW 2ds sI1 


where s = irpl = 7r/ 2 tan 2 9 e , and 


SLCJe g* n 

till Co \E\ 2 ds 111 


(6) 


( 7 ) 
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Fig. 5. Coordinate system for reflectarray efficiency analysis. 


where |/| 2 is the power distributed across the reflectarray aperture, II is the total power radiated by 
the feed within the reflectarray aperture, and III is the total power radiated by the feed to its front 
hemisphere. The analysis here assumes the feed has a cos' 7 # pattern and is polarized in the x direction. 
All the patch elements are also polarized in the x direction. The E vector in both Eqs. (6) and (7) is the 
feed far-field function and is given by 


E = 


\E\k = C 0 


cos 9 0 
r 



(8a) 


where the phase term e i kr is suppressed since phase does not play a role here for the efficiency calculation. 
The function P in Eq. (6) is the far-field pattern function of the patch element in the reflectarray and is 
assumed to be cos#. 


cos 9 6 


E * Px = C Q (cos # • cos 2 0 4- sin 2 0) cos 6 


(8b) 


The last term, cos#, in Eq. (8b) is the patch element pattern effect on the incoming wave from the feed. 

j _ f de f 2?r c ° ( cos<7+2 0 cos 2 0 4- cos^ 1 # sin 2 0) 

J o j O f 

Because ds = p dp d$,p = rsin0,r = // cos#, and dp = (f /cos 2 0)dd, 
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Because p = rsinO , 
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(ii) 


By using Eqs. (6), (7), (9), (10), and (11), 


Vilt = 


\I\ 2 _ [((1 ~ COS q + 1 <? e )/(7 + 1)) + ((1 - COS 9 g e )/<?)]' 

s77 “ 2 tan 2 6» e [(1 — cos 2 ^ 1 6> e )/(2qr + 1)] 


( 12 ) 


and 


Vs = JJJ = 1 - cos 29+1 9 e (13) 

By using Eqs. (12) and (13), the efficiencies are plotted against the feed pattern shapes in Fig. 6 with a 
given reflectarray diameter of 0.5 m at 32 GHz (Ka-band) and a f/D ratio of 1.0. The element spacing is 
half of the free-space wavelength, and the total number of patch elements is 8,937. These curves clearly 
indicate that the illumination and the spillover efficiencies are complementary to each other. The optimal 
aperture efficiency can be designed with a feed (] factor equal to 10.3. Another curve (Fig. 7) gives the 
optimum f/D ratio when an arbitrarily selected feed pattern is given (q = 8). Equations (12) and (13) 
are very important tools that aid in the design of an optimal reflectarray configuration. 

As mentioned previously, there are other efficiency factors, in addition to aperture efficiency, in the 
complete characterization of the reflectarray antenna system. The various contributors to the overall 
efficiency factor are estimated and listed in Table 1 for a Ka-band 0.5-m reflectarray with f/D - 1.0 and 
feed q factor = 10.3. 

It should be noted that the above feed loss does not include loss in the transmission line between the 
feed and the transceiver. 


V. Bandwidth Study 

Bandwidth is often an important quantity for satellite communication, especially with the increasing 
demand for higher data rates. At 32 GHz, 1 GHz of bandwidth (3 percent) is anticipated. In addition, 
it is expected that, for Ka-band deep space satellite communication, 32 GHz will be used for downlink 
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Fig. 6. Microstrip reflectarray spillover and illumination efficiencies versus 
feed pattern shape. 
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Table 1. Estimated efficiency of the microstrip refiectarray. 


Type of efficiency 

Efficiency, percent 

Loss, dB 

Illumination 

87 

0.60 

Spillover 

91 

0.41 

Patch loss 

97 

0.13 

Delay line loss 

95 

0.22 

Feed loss 

95 

0.22 

Cross-pol loss 

95 

0.22 

Total 

66 

1.80 


and 34 GHz for uplink. If a refiectarray antenna needs to cover both uplink and downlink frequencies, a 
bandwidth of more than 6 percent is required. Certainly, the bandwidth performance of the refiectarray 
is no match to the parabolic reflector, which has, theoretically, an infinite bandwidth. It is the intention 
of this section to study the bandwidth characteristics of the microstrip refiectarray and to optimize it for 
a given application. 

The bandwidth performance of a microstrip refiectarray can be limited by four factors: (1) the mi- 
crostrip patch element, (2) the array element spacing, (3) the feed antenna bandwidth, and (4) the 
differential spatial phase delay. Due to its thin cavity, the microstrip patch element can generally achieve 
a bandwidth of only 3 percent. To achieve a bandwidth larger than 3 percent, techniques such as the 
stacked dual patch or the patch with a thicker substrate can be employed. Ten- to fifteen- percent band- 
widths for microstrip antennas have been reported. The array element spacing limits the refiectarray 
performance such that, as frequency is decreased, the electrical element spacing becomes small, and ex- 
cessive mutual coupling effects start to degrade the array performance. As the frequency is increased, 
the electrical element spacing becomes large, and undesirable grating lobes begin to appear. Fortunately, 
previous calculations and experiences have shown that the element spacing effect will not be detrimental 
until the frequency variation is more than 30 percent (±15 percent around center frequency). The third 
bandwidth limiting factor is the feed antenna, which can be designed to operate over a bandwidth of at 
least 10 percent while maintaining a relatively constant beam shape and input impedance. Waveguide 
horns and cavity-backed dipoles are good candidates. If desired, an Archimedean spiral can be used to 
achieve more than 100 percent of bandwidth. The fourth limiting factor, differential spatial phase delay, 
has not been well understood and is separately detailed in the following paragraph. 

The differential spatial phase delay can best be explained by referring to Fig. 8, where the differential 
spatial phase delay, As, is the difference between the electrical paths 5i and S 2 . This As can be 
many multiples of the wavelength at the center operating frequency, such as As = NAX 0 , where N 
is an integer and A represents the fractional number of a free-space wavelength A 0 . At each patch 
location on the refiectarray, NA could be different numbers. In order to achieve constant aperture phase 
for the reradiated waves, the A at each patch location is compensated for by the appropriate length of the 
phase delay line attached to the patch. However, as frequency changes, the NA will change accordingly. 
Since the phase delay lines are fixed, a frequency excursion error will occur in the reradiated phase front. 
The old NAX 0 now becomes NA( X 0 + AA 0 ), where AA 0 is directly proportional to the frequency change. 
The amount of phase change is, therefore, NAAX 0 , which can be a significant portion of a wavelength 
(360 deg). To reduce the amount of frequency excursion error, the integer number N must be reduced. 
There are two ways to reduce N. One is to design the refiectarray with a larger f/D ratio, and the 
other is simply to use a refiectarray with a small electrical diameter. With a fixed f/D ratio, the larger 
the electrical diameter, the larger the N will be. The effects of the f/D ratio and the diameter on 
bandwidth performance are calculated by using Eq. (1) and are plotted in Figs. 9 and 10, where beam 
directivity versus frequency change is shown. The bandwidth effects of the patch element and the feed 
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Fig. 10. Antenna directivity versus frequency for a 1.0-m Ka-band microstrip reflectarray. 

antenna are not included in these figures. Figure 9 is plotted for a 32-GHz reflectarray having a diameter 
of 0.5 m and a total of 8,937 patch elements. Two f/D ratios of 0.5 and 1.0 are plotted in this figure. It 
is obvious that an f/D of 1.0 gives better bandwidth performance. Similar curves are plotted in Fig. 10 
for a 1-m-diameter reflectarray at 32 GHz. The number of patch elements in this case is 35,788. By 
comparing Figs. 9 and 10, the 1-m reflectarray with the same f/D gives less bandwidth than the 0.5-m 
one. This implies that, with the same f/D ratio, the larger the reflectarray, the smaller the bandwidth 
it will provide. The bandwidth performances of Figs. 9 and 10 are summarized in Table 2. 


Table 2. Bandwidth performance of the 
32-GHz reflectarray. 


f/D 

0.5-m diameter, 
percent 

1.0-m diameter, 
percent 

— 1-dB gain-drop bandwidth 

0.5 

4.8 

2.6 

1.0 

8.5 

4.5 

— 3-dB gain-drop bandwidth 

0.5 

8.4 

4.3 

1.0 

14.0 

7.5 
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The radiation patterns of the reflectarray will change as frequency changes. These differences are 
illustrated in Figs. 11 and 12 for two different f/D ratios. The pattern defocusing effect as frequency 
deviates from the center (design) frequency is clearly demonstrated in these figures. 


FYom the above results, it can be concluded that, among the four bandwidth limiting factors, the ele- 
ment spacing and the feed antenna are not serious concerns in designing the reflectarray if the bandwidth 
requirement is 15 percent or less. It also can be concluded that a 3-percent bandwidth (1-GHz bandwidth 
at 32 GHz) for the reflectarray is fairly easy to achieve. An 8-percent bandwidth to cover both the uplink 
frequency (34 GHz) and the downlink frequency (32 GHz) may require an f/D ratio close to 1.0 and a 
specially designed patch element. 


VI. Advantages of the Microstrip Reflectarray 

The microstrip reflectarray antenna takes the best characteristics and eliminates the poor features of 
the parabolic reflector, the array, and the microstrip patch. This contributes to the microstrip reflectar- 
ray’s six significant advantages, which are separately discussed below. 


A. Surface Mountable With Lower Mass and Volume 

Since the antenna’s reflecting surface is a thin, flat structure, it can be flush mounted onto the surface 
of a structure, such as a spacecraft’s main body or a building, with less supporting structure mass and 
volume as compared with the curved parabolic reflector. One possible application is shown in Fig. 13, 
where the reflectarray is surface mounted onto one side panel of a pentagonally shaped microspacecraft. 
Another application is given in Fig. 14, where the reflectarray is surface mounted onto a house wall for 
direct broadcast satellite (DBS) television service. In addition to the capability of surface mounting onto 
a flat structure, the reflectarray can also be mounted conformally onto a slightly curved structure (either 
concave or convex) . 1 he phase deviation of the curved structure can be compensated for in the design of 
the set of phase delay lines. 

B. Easily Deployable 

When a deployment mechanism is needed for a large aperture antenna, the flat structure of the 
reflectarray can be folded or unfolded by a simple hinge type of mechanism. A single or a double folding 
mechanism for a flat structure is significantly simpler (approximately an order of magnitude simpler) 
than that for a curved parabolic structure and is also more reliable. 1 The flat panel folding technique 
has been commonly used in the deployment of solar panels and has shown excellent reliability. 

C. Lower Manufacturing Cost 

The reflectarray, being in the form of a printed microstrip antenna, can be fabricated with a simple 
and low-cost etching process, especially when it is produced in large quantities. The antenna also can be 
cost effective just because of its flat structure. For example, the special molding process that is generally 
required for fabricating a curved paraboloid is not needed for a flat structure. 


D. Scannable Beam 

The main beam of the microstrip reflectarray can be designed to point at a large fixed angle (up to 
60 deg) from the broadside direction, while a parabolic reflector can only have a very limited beam tilt 
(several beamwidths). If cost permits, phase shifters can be placed in the phase delay transmission lines 
for electronic beam scanning. 


1 Personal communication with R. Freeland, Mechanical Engineer, Applied 
Laboratory, Pasadena, California, 1994. 
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E. Integratable With a Solar Array 

Since both the solar array and the reflectarray antenna are in the form of large, fiat panels, they have 
the possibility of being integrated together to save space and mass. There are two possible configurations 
for the integration. One is for a lower-frequency application (below 5 GHz) where the patch radiator is 
large enough to place solar cells on top of it. This is shown in Fig. 15, where the RF energy radiates 
from the perimeter of the patch and is not affected by the solar cells. The other configuration, proposed 
by Malibu Research [8], is to use thin wire dipoles as reflectarray radiators. This is illustrated in Fig. 16, 
where the thin dipoles are suspended on thin wires and placed in front of the solar panel. A metallic 
wire-mesh ground plane is placed between the dipoles and the solar cells. This ground plane is needed 
for the dipoles to reradiate effectively. This concept should be feasible for frequencies below 15 GHz. 
For frequencies higher than 15 GHz, the dipoles and the mesh ground plane may become too dense for 
significant amounts of sunlight to penetrate, and the mechanical support structure for the dipoles (with 
mesh ground plane) may be difficult to implement. Even for frequencies lower than 15 GHz, because the 
dipoles are placed in front of a wire mesh and not a continuous perfect ground plane, the RF efficiency 
will be impaired. In addition, because the dipoles and the wire mesh ground plane will partially block 
the sunlight, the solar cell efficiency will also be reduced. The degrees of impairment of RF and solar cell 
efficiencies remain to be studied. 

F. Very Large Aperture Array 

Due to the fact that no power divider is needed, the insertion loss of thousands of microstrip patches 
in the reflectarray is the same as that of a few patch elements. Thus, the reflectarray can achieve as good 
an efficiency as a large array antenna system. As a possible application, a very large, flat reflectarray 
can be constructed on flat land (array aperture parallel to the land) with its feed mounted on top of a 
high tower. The size of the antenna is only limited by the tower height. A feed height of 300 m with 
an f/D ratio of 0.5 will result in an array with a diameter of 600 m. At 10 GHz, this antenna size can 
produce a pencil beam 0.0035 deg in beamwidth and 95 dBi of directivity. Such an antenna can be used 
in astrophysics applications. 


VII. Proposed New Configurations of the Microstrip Reflectarray 

To improve the performance of the microstrip reflectarray, three novel concepts are introduced here. 
They are separately presented below. 

A. The Cassegrain Feed Reflectarray 

As mentioned previously, in order to achieve a wider bandwidth, a large f/D ratio is generally needed 
for the reflectarray. A large f/D implies that the focal feed has to protrude far from the array aperture, 
which will result in larger volume and larger mass. The proposed Cassegrain configuration, shown in 
Fig. 17, will reduce the feed height while maintaining the same or a higher effective f/D ratio. In 
addition, the transmission line loss between the feed and the transceiver is significantly reduced, which 
is especially important at higher frequencies, such as Ka-band. 

B. Bandwidth Enhancement by Subarray Arrangement Techniques 

The second concept uses specially arranged subarrays in the entire array aperture to improve the 
bandwidth performance. Figure 18 shows that, for linear polarization, each two-patch subarray uses the 
antiphase method [13] and, for circular polarization, each four-patch subarray uses the sequential rotation 
method [14] to improve the bandwidth. For conventional microstrip arrays, these subarraying techniques 
have been successfully demonstrated to increase bandwidth from 3 percent to more than 8 percent. For 
a single patch, as frequency changes away from its resonant frequency, higher order modes start to form 
in the patch’s cavity, and these will impair the input impedance match and polarization quality. By 
using the antiphase and sequential rotation methods, these high-order modes can be canceled, causing 
the desired fundamental mode to continue to dominate as the frequency changes. 
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Fig. 15. Integration of solar cells and microstrip radiators for low-frequency 
reflectarray applications. 


SOLAR PANEL 



Fig. 16. Integration of a solar panel with a thin-wire reflectarray antenna. 


C. Phase-Delay Compensation by Rotation of Identical Elements 

In a circularly polarized microstrip reflectarray, the circular polarization can be achieved by having 
two equal-length delay lines orthogonally attached to each patch (see Fig. 3). The length of these delay 
lines varies from element to element to compensate for the different spatial phase delays of the feed 
antenna. In this newly proposed technique, the required different phase-delay compensations for each 
of the elements are achieved by different angular rotations of the elements, as shown in Fig. 19. All 
elements, except for rotations, are identical with identical delay lines. The delay lines, in this case, only 
serve as the angular references for rotations and not for phase-delay compensations. Especially for a 
circular patch, without some kind of angular reference, the rotations between different patches cannot 
be differentiated. The required amount of rotation in degrees is half of that required for the phase delay 
in degrees. This technique of rotating the circularly polarized element to achieve the required phase has 
been demonstrated [15] previously in a conventional array. 



Fig. 18. Reflectarray bandwidth enhancement by subarray arrangement techniques: (a) linear 
polarization and (b) circular polarization. 


By combining the “rotation” technique with miniature motors, the main beam of the microstrip 
reflectarray can be made to scan. In other words, as shown in Fig. 20, a miniature motor can be 
placed under each patch element to mechanically rotate the patch to achieve the required phase for beam 
scanning. Since all the elements in a reflectarray are isolated from each other with no need for a power 
divider, no rotary joint is required, which will reduce cost and enhance reliability. By using the miniature 
motors, high-cost phase shifters are not needed for beam scanning. Consequently, a very large, phased 
array antenna system with relatively low cost may be realizable. 
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Fig. 19. Refiectarray phase-delay compensation by rotation of identical elements. 



Fig. 20. Mechanically phased refiectarray — beam scan by mechanical rotation of elements. 



VIII. Conclusion 

The circularly polarized microstrip reflectarray has been analyzed using conventional array theory. 
Antenna performance parameters, such as radiation pattern quality, directivity, efficiency, bandwidth, 
etc., are calculated for a 0.5-m Ka-band reflectarray with 8,937 patch elements. It is found that an 
8-percent bandwidth is achievable for this antenna and that bandwidth performance improves for a larger 
f/D ratio. An f/D ratio close to 1.0 is recommended for the 0.5-m Ka-band antenna. Overall antenna 
efficiency in the range of 50 to 70 percent is possible. Numerous advantages of this flat-plate reflector 
antenna have been discussed. Three novel configurations of the microstrip reflectarray are proposed for 
future studies. 
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Microwave performance testing of the new Deep Space Station (DSS)-24 34-m- 
diameter antenna was carried out during the summer of 1994. Efficiency measure- 
ments were made at the 8.45-GHz (X-band) and 32-GHz (Ka-band) frequencies 
both at the antenna Cassegrain (fl) and beam-waveguide (f3) focal points. In addi- 
tion, the antenna f3 efficiencies were measured on the DSS-24 operational 2.295-GHz 
(S-band) and 8.45-GHz feeds. This article presents the efficiency determinations as 
a function of elevation angle along with a corresponding error analysis of the mea- 
surements. Peak measured gains and efficiencies are tabulated for all frequencies. 


I. Introduction 

During an 8-week period in the summer of 1994, extensive microwave performance measurements 
were carried out on the new Deep Space Station (DSS)-24 34-m beam-waveguide (BWG) antenna at 
Goldstone, California. The testing period consisted of antenna efficiency and pointing calibration and 
system noise temperature measurements at the frequencies of 2.295 GHz (S-band), 8.45 GHz (X-band), 
and 32 GHz (Ka-band). 1 2 In addition, microwave holography and subsequent main antenna reflector-panel 
adjustments were carried out at the Cassegrain (fl) focus before the start of the efficiency calibration 
period. The X- and Ka-band efficiency measurements involved the use of portable microwave test packages 
installed at the fl focal point and the subterranean f3 beam- waveguide focal point. The test packages were 
previously developed for testing of the DSS-13 research and development (R.&D) antenna.' 1 In addition 
to the measurements on the R&D feeds, the performances of the DSS-24 operational S- and X-band feeds 
were also characterized at the f3 focal point. 

This article presents the efficiency calibration and analysis portion of the microwave DSS-24 testing. 
The methodology used to determine the antenna aperture efficiency from noise temperature measurements 
on radio source calibrators is reviewed and error analysis equations are presented. The following sections 
then present the results of the efficiency determination at the fl and f3 focal points, with discussion of 


1 The testing period is described in detail in DSS-24 Microwave Performance Characterization, Document 829-6 (internal 
document), Jet Propulsion Laboratory, Pasadena, California, May 15, 1994. 

2 M. J. Britcliffe, L. S. Alvarez, D. A. Bathker, P. W. Cramer, T. Y. Otoshi, D. J. Rochblatt, B. L. Seidel, S. D. Slobin, 
S. R. Stewart, W. Veruttipong, and G. E. Wood, DSS 13 Beam- Waveguide Antenna Project: Phase 1 Final Report , JPL 
D-8451 (internal document), Jet Propulsion Laboratory, Pasadena, California, May 15, 1991. 
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A disk-on-rod inside a corrugated horn is one of the horn configurations for 
dual-frequency or wide-band operation. A mode-matching analysis method is de- 
scribed. A disk-on-rod inside a corrugated horn is represented as a series of coaxial 
waveguide sections and circular waveguide sections connected to each other. Three 
kinds of junctions need to be considered: coaxial-to-coaxial, coaxial-to-circular, and 
circular-to-circular. A computer program was developed to calculate the scatter- 
ing matrix and the radiation pattern of a disk-on-rod inside a corrugated horn. 
The software was verified by experiment, and good agreement between calculation 
and measurement was obtained. The disk-on-rod inside a corrugated horn design 
gives an option to the Deep Space Network dual- frequency operation system , which 
currently is a two-horn/one-dichroic plate system. 


I. Introduction 

To design a dual-frequency horn for the DSS-13 beam waveguide antenna, an analysis tool needs 
to be developed. A side-view model of a circularly symmetric disk-on-rod inside a corrugated horn is 
shown in Fig. 1. The horn is subdivided to several sections that are either coaxial or circular wave- 
guide sections. The junctions between these sections are either coaxial-to-coaxial, coaxial-to-circular, or 
circular-to-circular waveguide junctions. In order to analyze the performance of a disk-on-rod inside a 
corrugated horn, a computer program based on the mode-matching method was developed [1,2]. The 
circular waveguide program and coaxial waveguide program that calculate the scattering matrix of the 
circular waveguide and coaxial waveguide of different sizes, respectively, are already available [3]. In order 
to simulate the disk-on-rod inside a corrugated horn, it was necessary to integrate the existing circular 
waveguide program and coaxial program with a new third program that handles the junction between 
circular and coaxial waveguides. 


II. Theory 

The analysis of the waveguide junctions is based on the mode-matching method. The following theory 
is used to calculate the scattering matrix of coaxial-to-circular waveguide junctions (Fig. 2). The elec- 
tromagnetic field is represented by coaxial waveguide modes in the coaxial waveguide region and circular 
waveguide modes in the circular waveguide region. Only waveguide modes of order 1 are considered in 
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pedestal room position. 9 The predicted antenna efficiency variation between the azimuth range of 60 to 
295 deg, computed for a perfectly aligned beam-waveguide mirror system, was 2 percent. As measured 
at DSS 24 and shown in Fig. 11, the variation is 4 percent over the same angle range. 


X- Summary of Results 

The aperture efficiency of the DSS-24 34-m beam-waveguide antenna has been calibrated at the S- 
(2.295-GHz), X- (8.45-GHz) and Ka-band (32-GHz) frequency bands. Table 2 presents a summary of the 
peak gains and efficiencies measured during the rf performance evaluation period. The uncertainties on 
the peak values are also given and are based on a propagation of measurement errors and radio-source 
flux and size correction errors. Table 3 gives the best-fit second-order polynomial coefficients for the X- 
and Ka-band efficiencies measured on the R&D feed packages. The first four peak gain and efficiency 
estimates in Table 2 were computed from these coefficients. The S- and X-band peak values measured 
on the operational feeds were derived by averaging the measurements in a small angular range about the 
observed peak efficiency angles. The Ka-band gravity loss as a function of elevation angle was computed 
from the best-fit DSS-24 f3 efficiency curve and is shown in Fig. 9. The final y- and z-axis subreflector 
(total) offset curves experimentally determined at the f3 focus are presented in Fig. 2 and have been left 
in the DSS-24 antenna-pointing controller for future tracking operations. 

Table 2. Summary of DSS-24 antenna peak gains and efficiencies at the S-, X-, 
and Ka-band frequencies. 


Frequency, GHz, 
and focal point 

Gain, dBi 

1-sigma gain 
error, dB 

Efficiency, 

percent 

1-sigma efficiency 
error, percent 

8.45 at fl 

68.33 

0.12 

75.25 

2.04 

8.45 at f3 

68.19 

0.11 

72.67 

1.81 

32 at fl 

78.97 

0.16 

60.60 

2.18 

32 at f3 

78.70 

0.20 

57.02 

2.51 

8.45 operational 

68.09 

0.14 

71.10 

2.30 

feed at f3 





2.295 operational 

56.79 

0.11 

71.50 

1.82 

feed at f3 a 






a Measured at azimuth angle = 180 deg. 


Table 3. Best-fit second-order polynomial coefficients for S- and Ka-band 
efficiencies, without atmosphere, measured on the R&D feed packages. 


Efficiency 

Efficiency model = co + c\el + C 2 eZ 2 
where el = elevation angle, deg 

Co 

c\ 

C2 

X-band at fl 

0.748067 

0.000211491 

-0.00000250421 

X-band at f3 

0.714972 

0.000457441 

-0.00000444267 

Ka-band at fl 

0.442223 

0.00737073 

-0.0000827938 

Ka-band at f3 

0.466983 

0.00504885 

-0.0000618488 


9 W. Veruttipong, “S-Band Efficiency and G/T Performance Performance for Various Azimuth Positions,” JPL Interoffice 
Memorandum 3327-94-113 (internal document), Jet Propulsion Laboratory, Pasadena, California, June 9, 1994. 
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The measurement errors, in general, are larger at f3 due to the decreased antenna gain and the 
apparent increase in the radio-source temperature measurement noise. The increased spread in the data 
in Figs. 6 and 7 is due to an azimuth dependence of the efficiency caused by a small beam-waveguide 
mirror misalignment. For the R&D feed location in the DSS-24 pedestal room (azimuth position = 
180 deg), the efficiency is slightly higher for sources with declinations greater than 35 deg rising and 
setting in the north. The curve fit in Fig. 6 is computed from the combined data sets of 3C274 and 3C84. 
The efficiency measurements on the northern-passing 3C84 (with peak elevation of 79 deg) almost all lie 
above the best-fit curve, while those of the southern-passing 3C274 (with peak elevation of 67 deg) are 
equal to or lie below it. 

A. Ka-Band Gravity Loss Analysis 

The gravity-induced roll-off of the f3 efficiency curve shown in Fig. 6 is flatter than that measured at 
fl (shown in Fig. 4). This is expected since the final subreflector offset curves determined at f3 with the 
automated subreflector optimization scheme provided better focusing than those used at fl. To quantify 
the performance increase, Fig. 8 presents the difference of the computed f3 and fl Ka-band gain curves 
after the f3 curve has been raised by 0.27 dBi to account for the difference in the measured peak gains. As 
seen, the low-elevation antenna gravity performance at f3 is superior to that at fl (e.g., 0.3 dBi greater at 
a 10-deg elevation). Thus, only the f3 efficiency curve should be used to predict the efficiency degradation 
as a function of elevation angle. The f3 gravity loss profile, computed from the best-fit efficiency curve, is 
shown in Fig. 9. It is assumed that there is zero loss at the peak gain elevation angle of 42 deg. Note that 
the asymmetry of the gravity loss curve is more probably due to an artifact of the differing profiles of the 
individual radio-source efficiency measurements, as discussed above, than to errors in the experimentally 
determined high-elevation angle subreflector focusing curves. 


VIII. Efficiency Measurements: X-Band Operational Feed 

Efficiency measurements at f3 were made on the X-band operational feed with the S-/X-band dichroic 
plate installed for the 3 days of 207-209. During this period, the computed zenith attenuation varied 
from 0.034 to 0.038 dB. The adjusted efficiency data points are shown in Fig. 10. 

On the X-band operational feed, the noise temperature measurements were taken at the input of the 
downconverter. However, as seen in Fig. 10, gain instability was a problem during the measurements. It 
was determined that the X-band maser was stable to ± 0.1 dB, corresponding to a 0.2-K noise temperature 
measurement stability on the principal calibrator, 3C274. Thus, at any elevation angle, the measurement 
noise on the efficiency measurements was 1.5 to 2 percent alone, rendering a quadratic curve fit impractical. 
The peak efficiency was estimated by averaging the data points between 42 and 48 deg. The peak efficiency 
on the X-band operational feed was computed to be 71.10 percent with an estimated uncertainty of 
2.30 percent. This corresponds to a peak X-band gain for the operational feed of 68.09 dBi ± 0.14 dBi. 


IX. Efficiency Measurements: S-Band Operational Feed 

S-band efficiency measurements were made on the operational feed on DOY 186-188. The radio sources 
tracked were 3C274 and 3C123, both southern-passing targets. Figure 11 shows the measurements, 
corrected for atmospheric absorption with an average zenith attenuation value of 0.03 dB, as a function 
of elevation angle. The azimuth angles for the data points range from 60 to 295 deg, with southern 
(180-deg azimuth) transit elevation angles of 67.5 deg for 3C274 and 83 deg for 3C123. The peak S- 
band efficiency was computed to be 71.50 percent, measured for both sources at the 180-deg azimuth 
angle. The uncertainty on the peak value is 1.82 percent, which corresponds to a peak S-band gain of 
56.79 dBi ± 0.11 dB. 

The variation of antenna efficiency as a function of azimuth angle shown in Fig. 11, with the peak 
operational gain at the 180-deg azimuth, was predicted by rf analysis for the S-band feed in the 270-deg 
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Fig. 5. DSS-24 32-GHz normalized efficiency uncertainty (a^ie/yri (el)) at the fl 
focal point, without atmosphere. 


temperature is maximum (i.e., the largest signal- to noise ratio condition). The uncertainty ratio increases 
when the elevation angle deviates from the rigging angle since the efficiency decreases due to gravity, and 
the largest magnitude is at the lower elevations due to the increased atmospheric attenuation variance 
a L(ei)- The one-standard deviation uncertainty <J n on the peak efficiency of 60.60 percent at 44.5 deg is 

2.18 percent (i.e., 60.6 x 0.036). The peak Ka-band efficiency and uncertainty at fl correspond to a peak 
gain of 78.97 dBi ± 0.16 dB. 


VII. Efficiency Measurements: Ka-Band at f3 

The Ka-band R&D microwave test package was moved to the f3 focal point after the X-band f3 
calibrations. The initial tracks focused on refining the subreflector offset curve determined at fl. The 
optimization procedure used at fl was automated in the automatic boresighting computer. Two days 
of star tracking yielded the final total subreflector-commanded position curves (labeled f3) shown in 
Fig. 2. The best-fit coefficients to the offset curves were input to an APC model file so that the offsets 
could be automatically applied during tracking operations. This is analogous to the pointing calibration 
and model determination. The final 2 days of tracking (DOY 174-175) provided the most stable and 
repeatable measurements for calibration. The zenith attenuation values for these 2 days, computed from 
SPC-10 weather recordings, varied from 0.10 to 0.12 dB, indicating a drier period than during the Ka- 
band measurements at fl. The best-fit efficiency curve and the adjusted data points for Ka-band at f3 
are shown in Fig. 6. 

The Ka-band peak efficiency at f3 was calculated to be 57.02 percent at an elevation angle of 42.0 deg. 
This is a decrease of 3.6 percent from the peak Ka-band fl efficiency. The uncertainty ratio, <r„ (e n/ij(ei), 
for the efficiency measurements is shown in Fig. 7. The uncertainty, <r v> on the peak efficiency of 57.02 per- 
cent at 42.0 deg is 2.51 percent. This corresponds to a peak gain of 78.70 dBi ± 0.20 dB. 
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EFFICIENCY COMPUTED ZENITH ATTENUATION VALUES, dB 



Fig. 3. Computed 32-GHz zenith attenuation values for DOY 158. 
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Fig. 4. DSS-24 32-GHz efficiency at the fl focal point, without atmosphere. 
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at most in an % 0.08-dB pointing loss for a 0.017-deg Ka-band full- width half-power beamwidth). A 
0.001-deg upper bound on the tracking error in each axis was a practical calibration goal given that 
the pointing-error measurement accuracy was not expected to be much better than 0.001 deg due to 
the low signal-to-noise operating conditions during radio source tracking. In addition to beam-pointing 
calibration, approximately 3 days were also devoted to obtaining the best subreflector offset positioning 
during the efficiency measurements. The optimization procedure, based on manually entering subreflector 
offsets into the APC in between boresights and postprocessing the output data together with the structural 
predict models, yielded the total fl offset curves shown in Fig. 2. These new focusing curves were then 
used for the remainder of the fl testing. The final 2 days of tracking (DOY 158 and 159) resulted in 
sufficient data points on 3C274 and 3C84 for the efficiency calibration curve. 
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Fig. 2. Subreflector position offset curves used during Ka-band 
measurements at fl and f3. 


The effects of the atmospheric attenuation become more significant with increasing frequency band. 
For this reason, the logged weather readings were closely reviewed. Figure 3 shows the computed 32-GHz 
zenith attenuation values for DOY 158, based on weather data logged at DSS 24 (sensors mounted on 
the X-band test package on the ground), and from readings logged daily at DSS 13 and SPC 10. In 
Fig. 3, the SPC-10 and DSS-13 values for A zen show good agreement, but the data points for DSS 24 
appear inconsistent, indicating a problem with the local sensors. Note that the readings at DSS 24 were 
logged only during the tracking periods. Since DSS 24 is closer in location to SPC 10 than to DSS 13, 
the SPC-10 weather data were used to remove the effects from the efficiency measurements. As seen in 
Fig. 3, the values for A zen varied from 0.11 to 0.16 dB during the course of the star tracks. The best-fit 
efficiency curve and the adjusted data points for Ka-band at fl are shown in Fig. 4. 

The Ka-band peak efficiency at fl was calculated to be 60.60 percent at an elevation angle of 44.5 deg. 
From Eq. (8), the uncertainty estimate for the efficiency measurements was computed for all elevation 
angles. A plot of the uncertainty ratio (^rj(ei)/v{ e ^) on efficiency data is shown in Fig. 5. As expected, 
the error in measuring the efficiency is minimum at 45-deg elevation, where the measured source noise 
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ELEVATION ANGLE, deg 

Fig. 1 . DSS-24 8.45-GHz efficiency at the fl and f3 focal points, without 

atmosphere. 

V, Efficiency Measurements: X-Band at f3 

X-band efficiency measurements at the f3 focal point were taken over 5 days (DOY 174-178) using the 
R&D microwave test package. Following the general data-taking strategy, the sources 3C274 and 3C84 
were tracked until efficiency curves could be computed on each source. The peak efficiency on the variable 
source 3C84 was then normalized to the peak 3C274 value. All of the measurements were compensated 
for atmospheric absorption on a point- by-point basis. For the final 2 days of tracking, the computed 
zenith attenuation values ranged only from 0.0334 to 0.035 dB. The best-fit efficiency curve based on 
these final days and the adjusted measurements is shown in Fig. 1. 

The peak X-band efficiency at f3 was calculated to be 72.67 percent at an elevation angle of 49 deg. 
From Eq. (8), the computed uncertainty on the peak efficiency is 1.81 percent. The peak X-band efficiency 
and uncertainty at f3 correspond to a peak gain of 68.19 dBi ±0.11 dB. As seen in Fig. 1, the efficiency 
data points taken on the two radio sources are numerous and repeatable over the elevation range. The 
uncertainty estimate of 1.81 percent is assumed constant with elevation angle since the last two error 
contributors (measurement noise and attenuation errors) in Eq. (8) contribute less than 10 percent to 
the total. As will be discussed below, radiometry techniques were applied to derive a total subreflector 
offset curve for use at Ka-band at f3 that differed from that used during the above X-band tracks and 
the Ka-band fl tracks. However, the effects on the presented X-band efficiency results are negligible. 

VI. Efficiency Measurements: Ka-Band at fl 

The Ka-band R&D microwave test package was used for efficiency calibration at both the DSS-24 fl 
and f3 focal points. At fl, just over a week was devoted to star tracks at Ka-band. The tracks from the 
first days, combined with the previous X-band boresight data, yielded the first pointing model for the 
antenna. This was requisite in order to keep the scan-to-scan error in each axis under 0.001 deg (resulting 
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III. Radio Source Calibration Values 

At all of the frequency bands, S-, X- and Ka-bands, the radio source 3C274 (Virgo A) is used as the 
principal calibrator. Due to a constrained time schedule during the radio frequency (rf) performance 
evaluation period, the main measurement strategy consisted of first obtaining a full track of 3C274 (rise 
to set) to establish the peak efficiency near the rigging angle (40-50 deg). Since 3C274 is only near a 
12-deg declination, higher declination sources (3C84 and 3C123) are then tracked to establish the upper 
efficiency curve. For the X- and Ka-bands, the efficiency curves generated for these latter sources are 
normalized to the measured 3C274 peak value. 

Table 1 summarizes the flux, source size correction values, calibration temperatures T a = 7i 0 o/C r , and 
the estimated uncertainties for the radio sources 3C274 and 30123 (the latter is used for calibration at 
S-band only). These values are taken directly from Richter 8 and are calculated for a DSN 34-m antenna 
with half-power beamwidths of 0.240 deg for S-band, 0.063 deg for X-band, and 0.017 deg for Ka-band. 
The uncertainties are used in the error analysis of the efficiency calculations. 


Table 1. Radio sources used for DSS-24 efficiency calibration. 


Source and 
frequency, GHz 

Flux density, 
Janskys 

C r 

Ta , K 

T a 1-sigma, K 

3C274 at 8.42 

44.69 

1.087 

13.518 

0.324 

3C274 at 32 

16.22 

1.270 

4.200 

0.130 

3C274 at 2.295 

137.24 

1.060 

42.570 

1.022 

3C123 at 2.295 

31.00 

1.000 

10.190 

0.245 


IV. Efficiency Measurements: X-Band at fl 

The X-band efficiency measurements were made at the fl focal point with the R&D feed package. As 
these were the first measurements made after holography, all of the star tracks applied the holography- 
determined best 45-deg subreflector offsets (in the y- and z-axis). These offsets are superimposed onto 
the structural (elevation-dependent) predict curves resident in the antenna-pointing controller (APC). No 
pointing model was used during the fl X-band tracks. However, subsequent analysis of the pointing errors 
(relative to the last boresight) indicated that no more than 0.002 deg of total pointing error was built up 
between the peak noise temperature estimates. This resulted in only a very small gain loss (~ 0.01 dB 
for a 0.063-deg X-band full-width half-power beamwidth) and, hence, measurements from the first two 
star tracks on 3C274 and 3C84 were used for efficiency calculation. 

The zenith attenuation values for the 2-day tracking period (DOY 144-145) varied from 0.036 to 
0.042 dB. The best- fit efficiency curve and the adjusted data points are shown in Fig. 1. The peak X- 
band efficiency at fl was determined to be 75.25 percent at an elevation angle of 42.5 deg. Using Eq. (8), 
the one-standard deviation uncertainty on the best-fit peak efficiency is computed to be 2.04 percent! 
The peak X-band efficiency and uncertainty correspond to a peak gain of 68.33 dBi ± 0.12 dB. The 
uncertainty on the efficiency below 30-deg elevation is greater due to insufficient (and nonrepeatable) 
measurement points available for the curve fit. At 20 deg, the one-standard deviation error is estimated 
to be 3 percent. 


8 P. Richter, op. cit. 
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T PiV (el) = L(el)T p (el) 


( 4 ) 


where T p (el) is the actual peak source temperature output from the boresight measurement, and the loss 
factor, L(el), is given by 


L{el) = 10 (>1(ei)/10) 


( 5 ) 


C. Error Analysis 

From Eq. (2), the variance of the efficiency, cr^, may be expressed as 



( 6 ) 


where or,, and a L are, respectively, the uncertainties of the absolute calibration source temperature, 
the measured source temperature, and the atmospheric attenuation loss factor. It can be shown that 
propagation of the loss factor, L(el), with respect to the zenith attenuation value, A zen , from Eq. (5) 
results in 


(^) 2 = ( 1 " ( 10 ) l ° s ( i ( e ,))) 2 fe ) 2 


which, inserted into Eq. (6), yields the efficiency variance function as 




+ (ln(10) log {L{el))f 



( 7 ) 


(8) 


A term-by-term computation of the error variance follows. The values of the pair (T a ,ar ft ) for various 
radio source calibrators are taken from Richter 6 for a 34-m antenna at each of the S-, X-, and Ka- 
frequency bands. The magnitude of the efficiency calibration is directly dependent on the value of r a 
and its accuracy. The uncertainty, <77^, is obtained from the residual sum of squares of the best-fit 
quadratic curve to the peak noise temperature samples, T p (e/), versus elevation angle. It is an estimate 
of the on-source radiometer measurement noise and is assumed to be constant with respect to time and 
elevation angle. The loss factor uncertainty, ctl( el) ) is computed from a probable error of 10 percent 
on the zenith attenuation value, 7 i.e., a Azen “ ®AA zen . Note that since A zen is calculated from surface 
weather measurements, it is a time- varying quantity. However, the error statistic <JA zen is modeled as 
stationary. For the DSS-24 calibration period, values for A zen are computed and analyzed using weather 
instrumentation recordings from DSS 24 (mounted on the microwave front-end test package itself), DSS 
13, and Signal Processing Center (SPC) 10. 

Lastly, it is assumed that there is no systematic efficiency loss due to beam-pointing errors during the 
boresight measurements since the algorithm is continually updating antenna-pointing corrections. Also, 
systematic efficiency losses due to errors in subreflector focusing and beam-waveguide mirror misalign- 
ments are, when applicable, added to the computed value of 


6 P. Richter, op. cit. 

7 Personal communication with S. Slobin, op. cit., August 1994. 
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the particular circumstances of each measurement set and the data analysis. The final section presents a 
summary of the peak gains and efficiencies measured during the testing period. 


II. Efficiency Calibration Methodology 

A. Efficiency Determination 

Determination of the DSS-24 aperture efficiency follows a previously developed methodology [l]. 3 Peak 
radio-source noise temperatures are measured at different antenna orientations via a boresight (actually a 
step-scan) algorithm [2] resident in a PC-based radiometer system. This automated boresighting program, 
know as AUTOBORE, was previously developed for performance testing of DSS 13 [3]. 

To yield an estimate of antenna efficiency, these measured quantities Eire then normalized by the 
absolute source calibration temperature, T a = T wo /C r , where Ti 00 is the 100-percent antenna efficiency 
temperature computed from the best-known radio source dux density and C r is the source size correction 
feictor (see [1] and [4] for more details). 4 

For this article, the aperture efficiency measured at DSS 24 is referenced to the input of the low-noise 
amplifier and, thus, includes the losses of the feed system. In equation form, the efficiency, T)(el ), for a 
given elevation angle, el, as measured on site is 


■q(el) = 


T P (el) 

T a 


( 1 ) 


where T p is the peak on-source temperature. The effects of atmospheric attenuation must be considered; 
then the Eq. (1) becomes 




(2) 


where L(el) is the atmospheric attenuation loss factor, which is described below. 

B. Correction for Atmospheric Attenuation 

In order to report the DSS-24 efficiency without the effects of atmospheric absorption, the estimated 
source peak temperatures from the antenna boresight measurements must be scaled up by a computed loss 
factor. A computer program, SDSATM4S.BAS, 5 is used to estimate total zenith attenuation values, A zen , 
from surface measurements of temperature, pressure, and relative humidity. From A zen , the attenuation 
loss factor for the logged boresight elevation angles can then be calculated. For clarity, the equations for 
the computation given in [1] are repeated here. The attenuation at elevation angle el is computed from 


A(el) = 


A 

** zen 

sin (el) 


(3) 


where A zen is in units of decibels. The radio source peak temperature, in a vacuum condition, is defined 

by 


3 Ibid. 

4 Values of the temperature T a = Tioo/CV for radio sources used to calibrate DSN antennas are given in P. Richter, 
DSN Radio Source List for Antenna Calibration , JPL D-3801, Rev. C (internal document), Jet Propulsion Laboratory,' 
Pasadena, California, August 19, 1993. 

5 Courtesy of S. Slobin, Telecommunications Systems Section, Jet Propulsion Laboratory, Pasadena, California. 
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CORRUGATED HORN 


DISK-ON-ROD . 



Fig. 1 . A model of a disk-on-rod inside a corrugated horn. 
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Fig. 2. Coaxial-to-circular waveguide junctions with (a) a different outer radius and 

(b) the same outer radius. 



the analysis, and the waveguide is _assumed to be nondissipative. In region I (the coaxial waveguide 
region), let the transverse field £7, at z = 0 be represented by the modal expansion 

M 

Ej = ^ {A m j T E m j) emi (1) 

m= 1 


M 

H i = ^ ^ {A m i — E m j ) h m j 
m= 1 


(2) 


where e m / and h m i are transverse modal field vectors and j4 m j,£ m / are the forward and reflected 
modal coefficients in the region I to be determined. The electromagnetic fields for the TE {e z = 0) and 
TM ( h z = 0) coaxial modes are as follows [4]. 

For the TE coaxial waveguide modes, the transverse e and h fields are 


e = e r f -f e<j)<p 


h = h r f + /i</>0 


(3) 

(4) 


er=c j_iMm sl „, 


(Xm r ) / b I 


00 — (jn £| 


,/ f X' m r\ 
l \b, ) 


(5) 


COS 0 


/ir — 




(6) 


L _ 

" r\TE 
' I rri 


where F\ is a combination of Bessel functions of the first kind, J v {z), and Bessel functions of the second 
kind, Y^z), of integral order v = 1, and is given by 


Fi 




(7) 


where Xm fh e root of the derivative of F\ when r = aj and C m is a normalization constant. The 
impedance of the TE waveguide mode, r ] TE , is given by 
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TE _ 


yV/g 

jyfvTe 

K VWAc) 2 -! 


for A < A c 
for A > A c 


( 8 ) 
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where n and £ are the absolute magnetic permeability and absolute dielectric constant of the medium, 
respectively. The free space wavelength is A, and A c is the cutoff wavelength (A c = A' m ) of the TE coaxial 
waveguide mode, which can be expressed as 


m = 1 

m = 2, 3,4, ■ • 

For the TM coaxial waveguide modes, the transverse e and h fields can be expressed as 

e r = - C m G[ cos<£ 


t 2 tt 

cm = «^) Tih^ {ai + bl) 
Km = «a//6/) - l)x^ (a/ “ b,) 


£<t> — On 


Cl {XmT/bl) 

X m r/b, 


sin (j> 


(9) 


( 10 ) 


( 11 ) 

h — er 
n<t> ~ „TM 
im 

where G\ is a combination of Bessel functions of the first kind, J v {x), and Bessel functions of the second 
kind, y„(x), of integral order v = 1, and is given by 

Gl (ir) “ J ' (tt) Y,fe ™ ) - Y (if) (12) 

where Xm is the mth root of G\ when r = aj. The impedance of the TM waveguide mode, r]™ , is 




TM 


\/?\A -(£) 2 for A < A c 

v r~ — ■ — 

~j\/^ ]/(-£) -1 for A > A c 


The cutoff wavelength of the TM coaxial waveguide mode, A c = A cm , can be expressed as 


(13) 


A 


cm — 


2tt 

(( aj/bj ) - 1 )x m 


(a/ - b[) 


(14) 


In region II, the transverse fields E n , Hu at 2 = 0 can be represented by the modal solution in region 
II (the circular waveguide region) as follows: 
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(15) 


N 

En = nil + B n u)e n n 

n= 1 

N 

Hji = y^(A n n - B n n)h n u (16) 

n= 1 


where e n i / and h n u are transverse modal fields and A n i[,B n jj are the forward and reflected modal 
coefficients in region II to be determined. 

The transverse e and h fields of the TE circular waveguide modes are 


■h ({Xn r )/ a Il) 
{Xn r )/ a H 


smcfr 


€$ — i 



h r = 


Vn E 



(17) 


(18) 


where \' n is the nth nonvanishing root of the derivative of the Bessel function J[(x’ n ) ~ 0- The C n is a 
normalization constant. The impedance of the TE waveguide mode, i } TE , is defined in Eq. (6), and the 
cutoff wavelength of the TE circular waveguide mode, A c = A(. n , can be expressed as 


v 2?r 

x cn = -y a n 
A.n 


(19) 


For TM circular waveguide modes, the transverse e and h fields are 


C r J 


&(j> — f-Ol 


Ji((x n r)/au) 


( 20 ) 
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and Xn is the nth nonvanishing root of Ji(Xn) = 0. The impedance of the TM waveguide mode, rj™ , 
is defined in Eq. (11), and the cutoff wavelength, A c = A cn , of the TM circular waveguide mode can be 
expressed as 


A 


27 r 

cn = &I I 

Xn 


( 22 ) 


By applying the boundary conditions that are discussed in detail in [1], the following pair of simultaneous 
matrix equations is obtained: 


+ [Bi]} = miAj,} + [B n ]} 


(23) 


{P] t {[Bu} -[>!//]} = [^{[A,] - [B,]} 


(24) 


where [,4/] and [B,\ are column matrices of M elements containing the unknown modal coefficients in 
region I; [>!//] and [ B n ] are column matrices of N elements containing the unknown modal coefficients 
in region II; [P} T is a transpose matrix of [P]; [P] is an M-by-N matrix; \Q] is an AT-by-JV diagonal 
matrix; and [R] is an M-by-M diagonal matrix. The elements of these three matrices, [PI, [Q] , and l R] 
are defined as follows: 


Qnn 


J Cml X IknII ’ ds 

S, 

J C-nll X hull ' ds 

Si i 


(25) 


(26) 
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mm 



x h mI • ds 


(27) 


In all cases, these integrals can be obtained in closed form. The Ftmm is the integration between two 
circular waveguide modes and Q nn is the integration between two coaxial waveguide modes in region I 
and region II, respectively; P mn is the integration between circular and coaxial waveguide modes and can 
be expressed as 
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(TM coaxial waveguide mode and TE circular waveguide mode) 

The submatrices [Sn], [S 12 ], [S 2 1 ], and [S 22 ] are derived from [P], [Q], and [P] by Eqs. (21) and (22): 


[Sn] = [VR] ([*] + [PflP ])- 1 ([«] - [P] T [^]) [VR] 


[s 12 ] = 2 [Vr] ([fl] + [P] T [P]) [P] T [Vq] 


-1 


[521] = 2 [v/q] (VQ + [ P ]{ P ] T ) [ P] T [^] “ 

[5 22 ] = [Vq] {[Q] + mT' ([Q] - [P}\P] T ) [Vq] 


(29) 

(30) 

(31) 

(32) 


When there are no dimension changes in the waveguide, which is then equivalent to a transmission 
line, the scattering matrix depends on the propagation constants of the waveguide modes in that straight 
section. To obtain the overall scattering matrix, the scattering matrices need to be cascaded. The 
procedure is described in detail in [1]. 

The above theory described the waveguide mode-matching method for coaxial-to-circular waveguide 
junctions. The same method was applied to circular-to-circular and coaxial-to-coaxial waveguide junctions 
(Figs. 3 and 4). 


III. Computer Program Development 

The geometrical configuration of the horn (flare angle, groove depth and width, and aperture size) 
is represented by a series of circular waveguides of different radii and lengths while the configuration of 
the disk-on-rod (diameter of the rod, diameter and thickness of the disk, disk spacing) is represented 
by a series of disks of various radii and lengths. These two data files are then combined and regenerated as 
a new data file containing the geometry of the disk-on-rod inside a horn. The new input data file in- 
cludes the outer radius (horn), inner radius (disk-on- rod), and the length of the (either coaxial or circular) 
waveguide section. The inner radius is zero for circular waveguide sections. All the sections are circu- 
larly symmetric with respect to the center axis of the feedhorn. The number of modes in each section is 
chosen according to the ratio M/N = (a/ - 6/)/(aj/ - &//) in order to converge to the correct values. 
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Fig. 3. Circular waveguide junctions. 
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Fig. 4. Coaxial waveguide junctions with a different (a) inner and outer 
radius, (b) outer radius, and (c) inner radius. 
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The P, Q, and R are calculated [Eqs. (25)— (27)] according to the type and size of the waveguide junctions. 
The scattering matrix of each section is calculated and cascaded with the scattering matrix of previous 
sections until the scattering matrix of the whole horn is obtained. 


The reflection matrix [Su] indicates the return loss of the horn. The radiation pattern can be computed 
by using the amplitude and phase of the transmitted circular waveguide modes ([S12]) at the horn aperture 
[5], Therefore, if the geometry of a feedhorn is available, both characteristics (radiation pattern and return 
loss) of the feedhorn are computed simply by inputting the geometrical dimensions to the computer 
program. 


IV. Verification of the Computer Program 

The program was first verified by comparing results with the existing circular waveguide program 
and coaxial waveguide program for appropriate junctions. Consistent results were achieved. Then an 
experiment was designed and performed in order to verify the scattering matrix of a circular-to-coaxial 
waveguide junction. The test piece is a WC137 circular waveguide with a circular aluminum rod suspended 
by two pieces of 0.0254 mm-thick kapton (Fig. 5). This structure includes a circular-to-coaxial and a 
coaxial-to-circular junction. The experiment was performed with rods of radii 7.62 and 10.16 mm and 
lengths of 63.5 and 76.2 mm. The amplitude and phase of the reflection (Su) and transmission (S21) 
coefficients were measured using a Hewlett Packard 8510C network analyzer. Good agreement between 
calculations and measurements was found in all the test cases. A test case of a rod of radius 7.62 mm 
and length 63.5 mm inside a WC137 circular waveguide of length 203.32 mm is shown in Figs. 6-9. The 
kapton in the waveguide was very thin in respect to the wavelength, so that it could be neglected in the 
computer modeling. 

The L-/C-band dual- frequency horn, which includes a C-band disk-on-rod inside a C-band launcher 
and an L-band horn, was also used to check the computer codes (Fig. 10) [6], By inputting the 
L-/C-band horn model in the program, the C-band and L-band radiation patterns were computed. 
Good agreement was shown between calculation and measurement at 5.01 GHz (C-band) and 1.668 GHz 
(L-band), respectively (Figs. 11-14) [7]. The slight asymmetry between the measured C-band E- and 
H-plane patterns was due to the C-band disk-on-rod becoming slightly off-centered during the trip to 
Goldstone, where the measurement was taken. All the results indicate that the software is reliable. 



Fig. S. The test piece for the circular-to-coaxial and coaxial-to-circular junction experiment. 
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Fig. 6. Calculated and measured amplitude of S^for a rod of radius 7.62 mm and length 63.5 mm 

inside a WC137 circular waveguide. 
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Fig. 7. Calculated and measured phase of S tJ for a rod of radius 7.62 mm and length 63.5 mm 

inside a WC137 circular waveguide. 
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Fig. 14. Measured and calculated E-plane pattern for the L-/C-band dual-frequency horn at 1.7 GHz. 
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V. Conclusion 

A disk-on-rod inside a horn was analyzed based on the mode-matching method. A computer program 
was developed to calculate the radiation pattern and the return loss of a horn by inputting the dimensions 
of the horn and the disk-on- rod. The computer program was verified by measurements and checked 
against other calculations. This software will be used to design an X-/Ka-band (8.45-GHz/33.7-GHz) 
dual-frequency horn for DSS 13. 
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